Análisis de Series Temporales

IOT Analytics

library(dplyr)
library(tidyverse)
library(kableExtra)
library(DT) # DataTables



library(imputeTS) # Tratamiento de datos faltantes

library(plotly)   # Gráficas
library(ggplot2)  # Gráficas
library(ggfortify)
library(gridExtra)

library(forecast) # Tratamiento de series temporales
load(file="Datos/DF_Energia_GMinutos.RData")
load(file="Datos/DF_Energia_GHoras.RData")
load(file="Datos/DF_Energia_GDiaria.RData")
load(file="Datos/DF_Energia_GMensual.RData")

Objetivo

Proyección del consumo energético a través de series temporales.

Procedimiento

Estudiaremos los datos de consumo energético correspondientes a los años 2007 a 2009 a través de series temporales.

Haremos proyección de futuro para el año 2010 y compararemos los resultados con los datos reales de este último año (serie real).

Series que vamos a estudiar:

  • Evolución mensual del consumo para los años 2007-2009. Predicción para el año 2010
  • Evolución diaria del consumo. Años 2007-2009. Predicción para el año 2010
  • Evolución semanal del consumo. Años 2007-2009 . Predicción para el año 2010

Evolución mensual del consumo. Años 2007-2009

Serie temporal y visualización

Vamos a representar la serie por meses para cada submedidor y tambien para la energía global co granularidad mensual. En este caso, la frecuencia será 12, tenemos una observación por mes para los años 2007,8 y 2009.

seriesMensuales <- Granularidad_meses %>% filter(year != 2006 & year != 2010)
## Creamos un objeto TS
tsSM1_Mensual<- ts(seriesMensuales$Sub_metering_1, frequency=12, start=c(2007,1))
tsSM2_Mensual<- ts(seriesMensuales$Sub_metering_2, frequency=12, start=c(2007,1))
tsSM3_Mensual<- ts(seriesMensuales$Sub_metering_3, frequency=12, start=c(2007,1))
tsSM_GAP_Mensual<- ts(seriesMensuales$Global_active_power, frequency=12, start=c(2007,1))
## Plot sub-meter 3 with autoplot - add labels, color
autoplot(tsSM1_Mensual, ts.colour = 'red', xlab = "Año", ylab = "Energía (Varios-hora)", main = "Evolución mensual del consumo de energía de la cocina \n Años 2007-2009")

autoplot(tsSM2_Mensual, xlab = "Año", ylab = "Energía (Varios-hora)", main = "Evolución mensual del consumo de energía de la lavandería \n Años 2007-2009")

autoplot(tsSM3_Mensual, ts.colour = 'blue', xlab = "Año", ylab = "Energía (Varios-hora)", main = "Evolución mensual del consumo de energía del aire acondicionado y termo eléctrico \n Años 2007-2009")

autoplot(tsSM_GAP_Mensual, ts.colour = 'green', xlab = "Año", ylab = "Energía (Varios-hora)", main = "Evolución mensual del consumo de energía \n Años 2007-2009")

Parece que los datos tienen estacionariedad (movimientos que se repiten cada año), esto ocurre para todas las variables.

Predicción del año 2010 antes de descomponer la serie

Vamos a plicar regresión lineal a cada serie temporal, y veremos el valor de \(R^2\) y del error cuadrático medio.

SM1. Cocina

Modelo

fitSM1_Mensual <- tslm(tsSM1_Mensual ~ trend + season) 
summary(fitSM1_Mensual)

Call:
tslm(formula = tsSM1_Mensual ~ trend + season)

Residuals:
     Min       1Q   Median       3Q      Max 
-20522.7  -5190.9    215.7   6726.8  16702.8 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  66531.9     7294.5   9.121 4.21e-09 ***
trend         -173.1      200.7  -0.862 0.397393    
season2     -19137.6     9635.6  -1.986 0.059061 .  
season3      -1263.5     9641.8  -0.131 0.896879    
season4     -15856.4     9652.3  -1.643 0.114034    
season5      -6134.4     9666.9  -0.635 0.531966    
season6      -9988.0     9685.6  -1.031 0.313160    
season7     -27046.6     9708.4  -2.786 0.010505 *  
season8     -38674.8     9735.4  -3.973 0.000602 ***
season9     -10641.7     9766.3  -1.090 0.287159    
season10    -15040.7     9801.3  -1.535 0.138536    
season11     -6105.3     9840.3  -0.620 0.541071    
season12     -3378.2     9883.2  -0.342 0.735595    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 11800 on 23 degrees of freedom
Multiple R-squared:  0.5796,    Adjusted R-squared:  0.3603 
F-statistic: 2.643 on 12 and 23 DF,  p-value: 0.02184

Resultados:

  • Valor del coeficiente de determinación: 0.5796. Es un valor bajo, el modelo únicamente explica el 58% de la variabilidad total.

  • Hay tres coeficientes significativamente no nulos.

Predicción IC

Predicción para los próximos 20 días

forecastfitSM1_Mensual <- forecast(fitSM1_Mensual, h=12, level=c(80,90))
## h=12 para predecir el año 2010
## Plot sub-meter 3 forecast, limit y and add labels
plot(forecastfitSM1_Mensual, ylab= "Watt-Hours", xlab="Fecha",
     main="Predicción del consumo energético del año 2010 \n  a través de un modelo de regresión")

SM2. Lavadero

Modelo

fitSM2_Mensual <- tslm(tsSM2_Mensual ~ trend + season) 
summary(fitSM2_Mensual)

Call:
tslm(formula = tsSM2_Mensual ~ trend + season)

Residuals:
   Min     1Q Median     3Q    Max 
-26325  -4033   1223   6470  16873 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  83389.3     7617.6  10.947 1.35e-10 ***
trend         -898.9      209.6  -4.289 0.000274 ***
season2     -13475.5    10062.5  -1.339 0.193598    
season3       7508.4    10069.0   0.746 0.463406    
season4     -15023.1    10079.9  -1.490 0.149708    
season5      -7727.5    10095.2  -0.765 0.451777    
season6     -11683.0    10114.7  -1.155 0.259931    
season7     -19075.8    10138.6  -1.882 0.072621 .  
season8     -28408.6    10166.7  -2.794 0.010305 *  
season9      -9204.4    10199.1  -0.902 0.376160    
season10      4809.1    10235.6   0.470 0.642890    
season11     -1696.0    10276.3  -0.165 0.870356    
season12     -3191.5    10321.1  -0.309 0.759941    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 12320 on 23 degrees of freedom
Multiple R-squared:  0.6541,    Adjusted R-squared:  0.4736 
F-statistic: 3.624 on 12 and 23 DF,  p-value: 0.003891

Resultados:

  • Valor del coeficiente de determinación: 0.6541. Es un valor bajo, el modelo únicamente explica el 65% de la variabilidad total.

  • Hay tres coeficientes significativamente no nulos.

Predicción IC

Predicción para los próximos 12 meses (año 2010)

forecastfitSM2_Mensual <- forecast(fitSM2_Mensual, h=12, level=c(80,90))
## h=12 para predecir el año 2010
## Plot sub-meter 3 forecast, limit y and add labels
plot(forecastfitSM2_Mensual, ylab= "Watt-Hours", xlab="Fecha",
     main="Predicción del consumo energético del año 2010 \n  a través de un modelo de regresión")

Se prevee una bajada del consumo

SM3. AC y termo

Modelo

fitSM3_Mensual <- tslm(tsSM3_Mensual ~ trend + season) 
summary(fitSM3_Mensual)

Call:
tslm(formula = tsSM3_Mensual ~ trend + season)

Residuals:
   Min     1Q Median     3Q    Max 
-87112 -18443   4781  17647  82107 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   299382      23189  12.911 5.06e-12 ***
trend           1866        638   2.926  0.00761 ** 
season2       -50805      30631  -1.659  0.11077    
season3       -29034      30651  -0.947  0.35336    
season4       -64190      30684  -2.092  0.04768 *  
season5       -54165      30730  -1.763  0.09125 .  
season6       -85197      30790  -2.767  0.01097 *  
season7      -144850      30863  -4.693 1.00e-04 ***
season8      -171144      30948  -5.530 1.27e-05 ***
season9       -69624      31047  -2.243  0.03486 *  
season10      -52696      31158  -1.691  0.10430    
season11      -37281      31282  -1.192  0.24550    
season12        7092      31418   0.226  0.82341    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 37510 on 23 degrees of freedom
Multiple R-squared:  0.7577,    Adjusted R-squared:  0.6312 
F-statistic: 5.992 on 12 and 23 DF,  p-value: 0.0001238

Resultados:

  • Valor del coeficiente de determinación: 0.7577. Es un valor que podría ser aceptable, el modelo únicamente explica el 76% de la variabilidad total.

  • Hay tres coeficientes significativamente no nulos.

Predicción IC

Predicción para los próximos 12 meses (año 2010)

forecastfitSM3_Mensual <- forecast(fitSM3_Mensual, h=12, level=c(80,90))
## h=12 para predecir el año 2010
## Plot sub-meter 3 forecast, limit y and add labels
plot(forecastfitSM3_Mensual, ylab= "Watt-Hours", xlab="Fecha",
     main="Predicción del consumo energético del año 2010 \n  a través de un modelo de regresión")

Se prevee una bajada del consumo

Comparación de los coeficientes y resultados de cada modelo

DatosRealesSeriesMensuales <- Granularidad_meses %>% filter(year == 2010 )

RMSE_SM1_GranMensual<-accuracy(forecastfitSM1_Mensual  ,DatosRealesSeriesMensuales$Sub_metering_1)[4]
R2_SM1_GranMensual<-0.5796
MASE_SM1_GranMensual<-accuracy(forecastfitSM1_Mensual  ,DatosRealesSeriesMensuales$Sub_metering_1)[12]

RMSE_SM2_GranMensual<-accuracy(forecastfitSM2_Mensual  ,DatosRealesSeriesMensuales$Sub_metering_2)[4]
R2_SM2_GranMensual<-0.6541
MASE_SM2_GranMensual<-accuracy(forecastfitSM2_Mensual  ,DatosRealesSeriesMensuales$Sub_metering_2)[12]



RMSE_SM3_GranMensual<-accuracy(forecastfitSM3_Mensual  ,DatosRealesSeriesMensuales$Sub_metering_3)[4]
R2_SM3_GranMensual<-0.7577
MASE_SM3_GranMensual<-accuracy(forecastfitSM3_Mensual  ,DatosRealesSeriesMensuales$Sub_metering_3)[12]

ResEvolMensual<-
cbind.data.frame(
      Submedidor=c("Cocina","Lavadero","AC y TermoE"),
      RMSE = c(RMSE_SM1_GranMensual,RMSE_SM2_GranMensual,RMSE_SM3_GranMensual),
     # MASE = c(MASE_SM1_GranMensual,MASE_SM2_GranMensual,MASE_SM3_GranMensual),
      R2 = c(R2_SM1_GranMensual,R2_SM2_GranMensual,R2_SM3_GranMensual)
)

ResEvolMensual %>% kable(booktabs=TRUE) %>% kable_styling(latex_options = "striped")
Submedidor RMSE R2
Cocina 11186.74 0.5796
Lavadero 10939.60 0.6541
AC y TermoE 57300.05 0.7577
# MASE: Mean absolute error
# RMSE: root mean scuare error. Diferencia entre los valores predichos y los observados.

El RMSE del modelo con mejor coeficiente de determinación es el más alto. Pero no hay que olvidar que el consumo eléctrico para el submedidor correspondiente al AC y termo eléctrico es el que más energía consume de todos con una diferencia muy grande.

Forecasting descomponiendo la serie (así ok)

Descomponer: quitar tendencia

Descomposición de la serie y visualización

Submetering 1

### Descomponer la serie
#Componentes_SM1_Mensual <- decompose(tsSM1_Mensual)
### Plot
#plot(Componentes_SM1_Mensual )
## stl(tsSM1_Mensual)
### Check summary statistics for decomposed sub-meter 1
#summary(Componentes_SM1_Mensual )

Se aprecia componente estacional: oscilacione que se repiten cada año. La componente tendencia tambien se aprecia. Tendencia decreciente.

 tsSM1_Mensual %>%
  stl( #t.window=13,
      s.window="periodic", 
      robust=TRUE)  %>% summary()
 Call:
 stl(x = ., s.window = "periodic", robust = TRUE)

 Time.series components:
    seasonal              trend            remainder        
 Min.   :-17333.003   Min.   :46768.55   Min.   :-37011.52  
 1st Qu.: -5016.432   1st Qu.:48940.13   1st Qu.: -3448.01  
 Median : -3280.804   Median :51834.57   Median :  -575.02  
 Mean   :     0.000   Mean   :51611.28   Mean   : -1053.45  
 3rd Qu.:  5593.844   3rd Qu.:53148.20   3rd Qu.:  3802.57  
 Max.   : 14381.497   Max.   :57632.43   Max.   : 32860.67  
 IQR:
     STL.seasonal STL.trend STL.remainder data 
     10610         4208      7251         16690
   %  63.6         25.2      43.4         100.0

 Weights:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.0000  0.9001  0.9452  0.8181  0.9965  1.0000 

 Other components: List of 5
 $ win  : Named num [1:3] 361 19 13
 $ deg  : Named int [1:3] 0 1 1
 $ jump : Named num [1:3] 37 2 2
 $ inner: int 1
 $ outer: int 15
tsSM1_Mensual %>%
  stl(t.window=13, s.window="periodic", robust=TRUE) %>%
  autoplot()

Remainder no es un muelle. Observamos tendencia creciente y fuerte componente estacional

Componentes_SM1_MensualSTL<- tsSM1_Mensual %>%stl( t.window=13,
  s.window="periodic",  robust=TRUE)
autoplot(Componentes_SM1_MensualSTL, ts.colour = 'blue')

plot.ts(tsSM1_Mensual, col = 'gray')
lines(Componentes_SM1_MensualSTL$time.series[,2], 
      col = "red", lwd = 1, lty = 2)

Submetering 2

## Descomponer la serie
Componentes_SM2_Mensual <- decompose(tsSM2_Mensual)
## Plot
plot(Componentes_SM2_Mensual )

## Check summary statistics for decomposed sub-meter 2
summary(Componentes_SM2_Mensual )
         Length Class  Mode     
x        36     ts     numeric  
seasonal 36     ts     numeric  
trend    36     ts     numeric  
random   36     ts     numeric  
figure   12     -none- numeric  
type      1     -none- character

Se aprecia componente estacional: oscilacione que se repiten cada año. La componente tendencia tambien se aprecia. Tendencia decreciente.

 tsSM2_Mensual %>%
  stl( #t.window=13,
      s.window="periodic", 
      robust=TRUE) %>% summary()
 Call:
 stl(x = ., s.window = "periodic", robust = TRUE)

 Time.series components:
    seasonal              trend            remainder        
 Min.   :-21630.641   Min.   :48453.75   Min.   :-43864.84  
 1st Qu.: -4761.576   1st Qu.:48533.62   1st Qu.: -5720.49  
 Median :  1499.661   Median :56646.72   Median :  -144.02  
 Mean   :    -0.001   Mean   :59759.83   Mean   : -1096.91  
 3rd Qu.:  5474.915   3rd Qu.:70783.48   3rd Qu.:  3351.28  
 Max.   : 15664.732   Max.   :77538.14   Max.   : 20892.53  
 IQR:
     STL.seasonal STL.trend STL.remainder data 
     10236        22250      9072         23422
   %  43.7         95.0      38.7         100.0

 Weights:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.0000  0.8448  0.9452  0.8851  0.9989  1.0000 

 Other components: List of 5
 $ win  : Named num [1:3] 361 19 13
 $ deg  : Named int [1:3] 0 1 1
 $ jump : Named num [1:3] 37 2 2
 $ inner: int 1
 $ outer: int 15
tsSM2_Mensual %>%
  stl(t.window=13, s.window="periodic", robust=TRUE) %>%
  autoplot()

Componentes_SM2_MensualSTL<- tsSM2_Mensual %>%stl( #t.window=13,
  s.window="periodic",  robust=TRUE)
autoplot(Componentes_SM2_MensualSTL, ts.colour = 'blue')

plot.ts(tsSM2_Mensual, col = 'gray')
lines(Componentes_SM2_MensualSTL$time.series[,2], 
      col = "red", lwd = 1, lty = 2)

Submetering 3

## ## Descomponer la serie
## Componentes_SM3_Mensual <- decompose(tsSM3_Mensual)
## ## Plot
## plot(Componentes_SM3_Mensual )
## ## Check summary statistics for decomposed sub-meter 2
## summary(Componentes_SM3_Mensual )

Se aprecia componente estacional: oscilacione que se repiten cada año. La componente tendencia tambien se aprecia. Tendencia creciente

tsSM3_Mensual %>%
  stl( #t.window=13,
      s.window="periodic", 
      robust=TRUE)  %>% summary()
 Call:
 stl(x = ., s.window = "periodic", robust = TRUE)

 Time.series components:
    seasonal              trend            remainder        
 Min.   :-107431.70   Min.   :249102.3   Min.   :-91715.37  
 1st Qu.: -12946.30   1st Qu.:250836.3   1st Qu.:-17758.23  
 Median :   9126.59   Median :276150.6   Median : -1364.69  
 Mean   :      0.00   Mean   :274156.0   Mean   : -2901.04  
 3rd Qu.:  25952.58   3rd Qu.:283886.9   3rd Qu.: 11936.12  
 Max.   :  68450.32   Max.   :316036.5   Max.   : 83277.64  
 IQR:
     STL.seasonal STL.trend STL.remainder data 
     38899        33051     29694         82030
   %  47.4         40.3      36.2         100.0

 Weights:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.0000  0.8028  0.9452  0.8268  0.9889  0.9999 

 Other components: List of 5
 $ win  : Named num [1:3] 361 19 13
 $ deg  : Named int [1:3] 0 1 1
 $ jump : Named num [1:3] 37 2 2
 $ inner: int 1
 $ outer: int 15
Componentes_SM3_MensualSTL<- tsSM3_Mensual %>%stl( #t.window=13,
  s.window="periodic",  robust=TRUE)
autoplot(Componentes_SM3_MensualSTL, ts.colour = 'blue')

plot.ts(tsSM3_Mensual, col = 'gray')
lines(Componentes_SM3_MensualSTL$time.series[,2], 
      col = "red", lwd = 1, lty = 2)

Predicción con Holt-Winters (suavizado exponencial) para el año 2010

Sin estacionalidad

Sería interesante verlo también con estacionalidad, para ver si se mantiene o no la tendencia.

Tendencia de predicción (crecimiento o decrecimiento del consumo) y predicción teniendo en cuenta las componentes

Evolución mensual del consumo (sin estacionalidad). Años 2007-2009

Cocina
#tsSM1_Mensual_Adjusted <- tsSM1_Mensual - Componentes_SM1_Mensual$seasonal
#autoplot(tsSM1_Mensual_Adjusted)
tsSM1_Mensual_detrended <- tsSM1_Mensual - Componentes_SM1_MensualSTL$time.series[,2]
plot.ts(tsSM1_Mensual_detrended, 
        ylab = "Energía (Vatios-hora)", 
        main = 'Seasonal variability of energy comsumption',
        cex.main = 0.85)

par(mfrow = c(1,2))
plot(Componentes_SM1_MensualSTL$time.series[,3], 
     col = "blue", main = 'Remainder', ylab = "")
qqnorm(Componentes_SM1_MensualSTL$time.series[,3])
qqline(Componentes_SM1_MensualSTL$time.series[,3])

shapiro.test(Componentes_SM1_MensualSTL$time.series[,3])

    Shapiro-Wilk normality test

data:  Componentes_SM1_MensualSTL$time.series[, 3]
W = 0.75719, p-value = 2.554e-06

Los resíduos no son solo ruido. Aún hay estacionalidad

Volvemos a descomponer la serie desestacionalizada una vez.

## Volvemos a descomponer la serie
plot(stl(tsSM1_Mensual_detrended ,s.window = "periodic",robust = TRUE))

La componente estacional se mueve en un ranfo de \(10^{-11}\), podemos asumir que se ha eliminado esta componente.

Suavizado exponencial de HoltWinters
tsSM1_HW_Mensual <- HoltWinters(tsSM1_Mensual_detrended ,
## Holt Winters Exponential Smoothing & Plot
                              #  tsSM1_Mensual_Adjusted, 
  beta=TRUE, gamma=TRUE)

# optimiza
plot(tsSM1_HW_Mensual,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Pronóstico de HoltWinters para el año 2010

Vamos a predecir las próximas 12 observaciones.

## HoltWinters forecast & plot
tsSM1_HW_Mensual_for <- forecast(tsSM1_HW_Mensual, h=12)
plot(tsSM1_HW_Mensual_for , ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico",
     main="Predicción del consumo energético de la cocina \ Año 2010" ) 

## Forecast HoltWinters with diminished confidence levels
tsSM1_HW_Mensual_forC <- forecast(tsSM1_HW_Mensual, h=12, level=c(10,25))
## Plot only the forecasted area
plot(tsSM1_HW_Mensual_forC, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico", main="Predicción para el año 2010", start(2010))

Lavadero
## Seasonal adjusting sub-meter 3 by subtracting the seasonal component & plot
# tsSM2_Mensual_Adjusted <- tsSM2_Mensual - Componentes_SM2_Mensual$seasonal
# autoplot(tsSM2_Mensual_Adjusted)
tsSM2_Mensual_detrended <- tsSM2_Mensual - Componentes_SM2_MensualSTL$time.series[,2]
plot.ts(tsSM2_Mensual_detrended, 
        ylab = "Energía (Vatios-hora)", 
        main = 'Seasonal variability of energy comsumption',
        cex.main = 0.85)

par(mfrow = c(1,2))
plot(Componentes_SM2_MensualSTL$time.series[,3], 
     col = "blue", main = 'Remainder', ylab = "")
qqnorm(Componentes_SM2_MensualSTL$time.series[,3])
qqline(Componentes_SM2_MensualSTL$time.series[,3])

shapiro.test(Componentes_SM2_MensualSTL$time.series[,3])

    Shapiro-Wilk normality test

data:  Componentes_SM2_MensualSTL$time.series[, 3]
W = 0.85783, p-value = 0.0002898

Los resíduos no siguen una distribución normal. La descomposición no es buena.

## Volvemos a descomponer la serie
plot(stl(tsSM2_Mensual_detrended ,s.window = "periodic",robust = TRUE))

La componente estacional se mueve en un ranfo de \(10^{-11}\), podemos asumir que se ha eliminado esta componente.

Suavizado exponencial de HoltWinters
## Holt Winters Exponential Smoothing & Plot
tsSM2_HW_Mensual <- HoltWinters(tsSM2_Mensual_detrended ,
                                # tsSM2_Mensual_Adjusted,
                                beta=TRUE, gamma=TRUE)
plot(tsSM2_HW_Mensual,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Pronóstico de HoltWinters para el año 2010

Vamos a predecir las próximas 12 observaciones.

## HoltWinters forecast & plot
tsSM2_HW_Mensual_for <- forecast(tsSM2_HW_Mensual, h=12)
plot(tsSM2_HW_Mensual_for , ylab= "Vatios-Hora", xlab="Fecha - Lavadero",
     main="Predicción del consumo energético del lavadero \ Año 2010" ) 

## Forecast HoltWinters with diminished confidence levels
tsSM2_HW_Mensual_forC <- forecast(tsSM2_HW_Mensual, h=12, level=c(10,25))
## Plot only the forecasted area
plot(tsSM2_HW_Mensual_forC, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico", main="Predicción para el año 2010", start(2010))

AC y termo eléctrico
# ## Seasonal adjusting sub-meter 3 by subtracting the seasonal component & plot
# tsSM3_Mensual_Adjusted <- tsSM3_Mensual - Componentes_SM3_Mensual$seasonal
# autoplot(tsSM3_Mensual_Adjusted)
tsSM3_Mensual_detrended <- tsSM3_Mensual - Componentes_SM3_MensualSTL$time.series[,2]
plot.ts(tsSM3_Mensual_detrended, 
        ylab = "Energía (Vatios-hora)", 
        main = 'Seasonal variability of energy comsumption',
        cex.main = 0.85)

par(mfrow = c(1,2))
plot(Componentes_SM3_MensualSTL$time.series[,3], 
     col = "blue", main = 'Remainder', ylab = "")
qqnorm(Componentes_SM3_MensualSTL$time.series[,3])
qqline(Componentes_SM3_MensualSTL$time.series[,3])

shapiro.test(Componentes_SM3_MensualSTL$time.series[,3])

    Shapiro-Wilk normality test

data:  Componentes_SM3_MensualSTL$time.series[, 3]
W = 0.93643, p-value = 0.03938

Los resíduos no siguen una distribución normal. La descomposición no es buena.

## Volvemos a descomponer la serie
plot(stl(tsSM3_Mensual_detrended ,s.window = "periodic",robust = TRUE))

La componente estacional se mueve en un ranfo de \(10^{-11}\), podemos asumir que se ha eliminado esta componente.

Suavizado exponencial de HoltWinters
## Holt Winters Exponential Smoothing & Plot
tsSM3_HW_Mensual <- HoltWinters(tsSM3_Mensual_detrended ,
                              # tsSM3_Mensual_Adjusted,
                                beta=TRUE, gamma=TRUE)
plot(tsSM3_HW_Mensual,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Pronóstico de HoltWinters para el año 2010

Vamos a predecir las próximas 12 observaciones.

## HoltWinters forecast & plot
tsSM3_HW_Mensual_for <- forecast(tsSM3_HW_Mensual, h=12)
plot(tsSM3_HW_Mensual_for , ylab= "Vatios-Hora", xlab="Fecha - Lavadero",
     main="Predicción del consumo energético del AC y termo eléctrico \ Año 2010" ) 

## Forecast HoltWinters with diminished confidence levels
tsSM3_HW_Mensual_forC <- forecast(tsSM3_HW_Mensual, h=12, level=c(10,25))
## Plot only the forecasted area
plot(tsSM3_HW_Mensual_forC, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico", main="Predicción para el año 2010", start(2010))

Comparación de los coeficientes y resultados de cada predicción. Datos sin estacionalidad

En las gráficas de predicción, hemos observado que el ajuste no mantiene la tendencia de consumo energético de la serie.

RMSE_SM1_GranMensualFor<-accuracy(tsSM1_HW_Mensual_forC   ,DatosRealesSeriesMensuales$Sub_metering_1)[4]
RMSE_SM2_GranMensualFor<-accuracy(tsSM2_HW_Mensual_forC   ,DatosRealesSeriesMensuales$Sub_metering_2)[4]
RMSE_SM3_GranMensualFor<-accuracy(tsSM3_HW_Mensual_forC   ,DatosRealesSeriesMensuales$Sub_metering_3)[4]


ResEvolMensual2_SinEstacionalidadConTendencia <-
cbind.data.frame(
      Submedidor=c("Cocina","Lavadero","AC y TermoE"),
      RMSE = c(RMSE_SM1_GranMensualFor,
               RMSE_SM2_GranMensualFor,
               RMSE_SM3_GranMensualFor)
)
ResEvolMensual2_SinEstacionalidadConTendencia %>% kable(booktabs=TRUE) %>% kable_styling(latex_options = "striped")
Submedidor RMSE
Cocina 45122.42
Lavadero 44294.27
AC y TermoE 318158.49

Evolución mensual del consumo (con estacionalidad). Años 2007-2009

Vamos a ver una predicción del consumo enegético de cada submedidor con el método de HW. En este caso, no eliminaremos la estacionalidad, para ver así si se mantiene la tendencia de consumo.

Cocina
tsSM1_Mensual_AdjustedII <- tsSM1_Mensual # - Componentes_SM1_Mensual$seasonal
autoplot(tsSM1_Mensual_AdjustedII)

Parece un muelle.

## Volvemos a descomponer la serie
plot(decompose(tsSM1_Mensual_AdjustedII))

La componente estacional se mueve en un rango de \(10^{-11}\), podemos asumir que se ha eliminado esta componente.

Suavizado exponencial de HoltWinters
## Holt Winters Exponential Smoothing & Plot
tsSM1_HW_MensualII <- HoltWinters(tsSM1_Mensual_AdjustedII, beta=TRUE, gamma=TRUE)
# optimiza
plot(tsSM1_HW_MensualII,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Pronóstico de HoltWinters para el año 2010

Vamos a predecir las próximas 12 observaciones.

## HoltWinters forecast & plot
tsSM1_HW_Mensual_forII <- forecast(tsSM1_HW_MensualII, h=12)
plot(tsSM1_HW_Mensual_forII , ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico",
     main="Predicción del consumo energético de la cocina \ Año 2010" ) 

## Forecast HoltWinters with diminished confidence levels
tsSM1_HW_Mensual_forCII <- forecast(tsSM1_HW_MensualII, h=12, level=c(10,25))
## Plot only the forecasted area
plot(tsSM1_HW_Mensual_forCII, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico", main="Predicción para el año 2010", start(2010))

Lavadero
## Seasonal adjusting sub-meter 3 by subtracting the seasonal component & plot
tsSM2_Mensual_AdjustedII <- tsSM2_Mensual ## - Componentes_SM2_Mensual$seasonal
autoplot(tsSM2_Mensual_AdjustedII)

No arece un muelle.

## Volvemos a descomponer la serie
plot(decompose(tsSM2_Mensual_AdjustedII)) 

TENDENCIA CLARAMENTE DECRECIENTE. Tambien observamos esacionalidad

Suavizado exponencial de HoltWinters
## Holt Winters Exponential Smoothing & Plot
tsSM2_HW_MensualII <- HoltWinters(tsSM2_Mensual_AdjustedII, beta=TRUE, gamma=TRUE)
plot(tsSM2_HW_MensualII,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Pronóstico de HoltWinters para el año 2010

Vamos a predecir las próximas 12 observaciones.

## HoltWinters forecast & plot
tsSM2_HW_Mensual_forII <- forecast(tsSM2_HW_MensualII, h=12)
plot(tsSM2_HW_Mensual_forII , ylab= "Vatios-Hora", xlab="Fecha - Lavadero",
     main="Predicción del consumo energético del lavadero \ Año 2010" ) 

## Forecast HoltWinters with diminished confidence levels
tsSM2_HW_Mensual_forCII <- forecast(tsSM2_HW_MensualII, h=12, level=c(10,25))
## Plot only the forecasted area
plot(tsSM2_HW_Mensual_forCII, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico", main="Predicción para el año 2010", start(2010))

AC y termo eléctrico
## Seasonal adjusting sub-meter 3 by subtracting the seasonal component & plot
tsSM3_Mensual_AdjustedII <- tsSM3_Mensual # - Componentes_SM3_Mensual$seasonal
autoplot(tsSM3_Mensual_AdjustedII)

No parece un muelle.

## Volvemos a descomponer la serie
plot(decompose(tsSM3_Mensual_AdjustedII))

Suavizado exponencial de HoltWinters
## Holt Winters Exponential Smoothing & Plot
tsSM3_HW_MensualII <- HoltWinters(tsSM3_Mensual_AdjustedII, beta=TRUE, gamma=TRUE)
plot(tsSM3_HW_MensualII,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Pronóstico de HoltWinters para el año 2010

Vamos a predecir las próximas 12 observaciones.

## HoltWinters forecast & plot
tsSM3_HW_Mensual_forII <- forecast(tsSM3_HW_MensualII, h=12)
plot(tsSM3_HW_Mensual_forII , ylab= "Vatios-Hora", xlab="Fecha - Lavadero",
     main="Predicción del consumo energético del AC y termo eléctrico \ Año 2010" ) 

## Forecast HoltWinters with diminished confidence levels
tsSM3_HW_Mensual_forCII <- forecast(tsSM3_HW_MensualII, h=12, level=c(10,25))
## Plot only the forecasted area
plot(tsSM3_HW_Mensual_forCII, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico", main="Predicción para el año 2010", start(2010))

Comparación de los coeficientes y resultados de cada predicción. Datos con estacionalidad
RMSE_SM1_GranMensualForII<-accuracy(tsSM1_HW_Mensual_forCII   ,DatosRealesSeriesMensuales$Sub_metering_1)[4]
RMSE_SM2_GranMensualForII<-accuracy(tsSM2_HW_Mensual_forCII   ,DatosRealesSeriesMensuales$Sub_metering_2)[4]
RMSE_SM3_GranMensualForII<-accuracy(tsSM3_HW_Mensual_forCII   ,DatosRealesSeriesMensuales$Sub_metering_3)[4]


ResEvolMensual2SinTendenciaConEstacionalidad<-
cbind.data.frame(
      Submedidor=c("Cocina","Lavadero","AC y TermoE"),
      RMSE = c(RMSE_SM1_GranMensualForII,
               RMSE_SM2_GranMensualForII,
               RMSE_SM3_GranMensualForII)
   
    
)

ResEvolMensual2SinTendenciaConEstacionalidad %>% kable(booktabs=TRUE) %>% kable_styling(latex_options = "striped")
Submedidor RMSE
Cocina 16362.16
Lavadero 27223.90
AC y TermoE 165249.73

Comparación de resultados. Predicción antes y despues de eliminar estacionalidad

Comparacion<-cbind.data.frame(
  ResEvolMensual2_SinEstacionalidadConTendencia,
  ResEvolMensual2SinTendenciaConEstacionalidad[,2])

colnames(Comparacion) <- c("Submedidor","RMSE Con Estacionalidad","RMSE Sin estacionalidad")
Comparacion %>% kable(booktabs=TRUE) %>% kable_styling(latex_options = "striped")
Submedidor RMSE Con Estacionalidad RMSE Sin estacionalidad
Cocina 45122.42 16362.16
Lavadero 44294.27 27223.90
AC y TermoE 318158.49 165249.73

Resultados más fiables y con menor error trás haber eliminado la estacionalidad de la serie. También así podemos ver la tendencia.

#p1<-autoplot(tsSM1_Mensual_Adjusted)
#p2<-autoplot(tsSM1_Mensual_AdjustedII)
#grid.arrange(
#             grobs = list(p1,p2),
#           nrow = 2,
#             top="Sin estacionalidad",bottom="Con estacionalidad"
#             )

Comparación de resultados. Predicción antes y despues de eliminar la tendencia (descomponer)

En las gráficas de predicción de suaviazado exponencial, hemos observado que el ajuste cuando eliminamos la estacionalidad, no mantiene la tendencia de consumo energético de la serie. Sin embargo, al mantener la estacionalidad se puede apreciar como las predicciones ajustan la tendencia de la serie original

# Predicción antes y despues de eliminar estacionalidad:
Comparacion
   Submedidor RMSE Con Estacionalidad RMSE Sin estacionalidad
1      Cocina                45122.42                16362.16
2    Lavadero                44294.27                27223.90
3 AC y TermoE               318158.49               165249.73
# Sin descomponer


ComparacionII<-cbind.data.frame(ResEvolMensual,
                              Comparacion[-1]
                              )
colnames(ComparacionII) <- c("Submedidor","RMSE","R^2","Con estacionalidad","Sin estacionalidad")
ComparacionII %>% kable(booktabs=TRUE) %>% 
  kable_styling(latex_options = "striped") %>% 
  add_header_above( c("", "Serie sin descomponer" = 2, "Descomponiendo la serie" = 2)) %>% 
  add_header_above( c("", "RMSE" = 4))
RMSE
Serie sin descomponer
Descomponiendo la serie
Submedidor RMSE R^2 Con estacionalidad Sin estacionalidad
Cocina 11186.74 0.5796 45122.42 16362.16
Lavadero 10939.60 0.6541 44294.27 27223.90
AC y TermoE 57300.05 0.7577 318158.49 165249.73

Evolución semanal del consumo. Años 2007-2009

Serie temporal y visualización

Vamos a representar la serie por meses para cada submedidor y tambien para la energía global co granularidad mensual. En este caso, la frecuencia será 12, tenemos una observación por mes para los años 2007,8 y 2009.

Granularidad_semanas <- Granularidad_horas %>% group_by(year, weekday) %>% 
   summarize(
      Sub_metering_1 = sum (Sub_metering_1),
      Sub_metering_2 = sum (Sub_metering_2),
      Sub_metering_3 = sum (Sub_metering_3),
      Global_active_power = sum (Global_active_power),
      energia2 = sum (energia2)
   )
seriesSemanales <- rbind.data.frame(
   filter(Granularidad_semanas, year== 2007  ),
   filter(Granularidad_semanas, year== 2008  ),
   filter(Granularidad_semanas, year== 2009  )
)

## Creamos un objeto TS
   tsSM1_Semanal<- ts(seriesSemanales$Sub_metering_1,         frequency=7, start=c(2007,1))
   tsSM2_Semanal<- ts(seriesSemanales$Sub_metering_2,         frequency=7, start=c(2007,1))
   tsSM3_Semanal<- ts(seriesSemanales$Sub_metering_3,         frequency=7, start=c(2007,1))
tsSM_GAP_Semanal<- ts(seriesSemanales$Global_active_power, frequency=7, start=c(2007,1))
## Plot sub-meter 3 with autoplot - add labels, color
autoplot(tsSM1_Semanal, ts.colour = 'red', xlab = "Semana del año", ylab = "Energía (Varios-hora)", main = "Evolución semanal del consumo de energía de la cocina \n Años 2007-2009")

autoplot(tsSM2_Semanal, xlab = "Semana del año", ylab = "Energía (Varios-hora)", main = "Evolución semanal del consumo de energía de la lavandería \n Años 2007-2009")

autoplot(tsSM3_Semanal, ts.colour = 'blue', xlab = "Semana del año", ylab = "Energía (Varios-hora)", main = "Evolución semanal del consumo de energía del aire acondicionado y termo eléctrico \n Años 2007-2009")

autoplot(tsSM_GAP_Semanal, ts.colour = 'green', xlab = "Semana del año", ylab = "Energía (Varios-hora)", main = "Evolución semanal del consumo de energía \n Años 2007-2009")

AC y Termo: podría apreciarse cierta estacionalidad Energía total: muelle

Predicción del año 2010 antes de descomponer la serie

Vamos a plicar regresión lineal a cada serie temporal, y veremos el valor de \(R^2\) y del error cuadrático medio.

SM1. Cocina

Modelo

fitSM1_Semanal <- tslm(tsSM1_Semanal ~ trend + season) 
summary(fitSM1_Semanal)

Call:
tslm(formula = tsSM1_Semanal ~ trend + season)

Residuals:
   Min     1Q Median     3Q    Max 
-11708  -2978   1845   4714   7123 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  59656.6     4512.2  13.221 6.48e-09 ***
trend         -508.6      263.8  -1.928  0.07593 .  
season2      12293.0     5646.7   2.177  0.04850 *  
season3      23297.9     5665.1   4.113  0.00122 ** 
season4       7005.5     5695.7   1.230  0.24051    
season5      13062.8     5738.3   2.276  0.04039 *  
season6      81428.1     5792.6  14.057 3.07e-09 ***
season7      91174.1     5858.3  15.563 8.75e-10 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 6908 on 13 degrees of freedom
Multiple R-squared:  0.975, Adjusted R-squared:  0.9616 
F-statistic: 72.49 on 7 and 13 DF,  p-value: 2.122e-09

Resultados:

  • Valor del coeficiente de determinación: 0.975. Es un valor bajo, el modelo únicamente explica el 97.5% de la variabilidad total.

  • Hay tres coeficientes significativamente no nulos.

Predicción IC

Predicción para los próximos 20 días

forecastfitSM1_Semanal <- forecast(fitSM1_Semanal, h=7, level=c(80,90))
## h=12 para predecir el año 2010
## Plot sub-meter 3 forecast, limit y and add labels
plot(forecastfitSM1_Semanal, ylab= "Watt-Hours", xlab="Fecha",
     main="Predicción del consumo energético semanal total del año 2010 \n  a través de un modelo de regresión")

SM2. Lavadero

Modelo

fitSM2_Semanal <- tslm(tsSM2_Semanal ~ trend + season) 
summary(fitSM2_Semanal)

Call:
tslm(formula = tsSM2_Semanal ~ trend + season)

Residuals:
   Min     1Q Median     3Q    Max 
-34037 -13452  -3470   8762  38539 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  88194.9    15111.4   5.836 5.82e-05 ***
trend        -2641.6      883.3  -2.991 0.010427 *  
season2      43786.9    18910.7   2.315 0.037565 *  
season3      73582.1    18972.5   3.878 0.001903 ** 
season4      -4120.0    19075.0  -0.216 0.832350    
season5      28106.0    19217.7   1.463 0.167355    
season6      64649.2    19399.5   3.333 0.005399 ** 
season7      83987.4    19619.5   4.281 0.000895 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 23140 on 13 degrees of freedom
Multiple R-squared:  0.7717,    Adjusted R-squared:  0.6487 
F-statistic: 6.276 on 7 and 13 DF,  p-value: 0.002282

Resultados:

  • Valor del coeficiente de determinación: 0.7717. Es un valor aceptable, el modelo únicamente explica el 77.17% de la variabilidad total.

  • Hay tres coeficientes significativamente no nulos.

Predicción IC

Predicción para los próximos 12 meses (año 2010)

forecastfitSM2_Semanal <- forecast(fitSM2_Semanal, h=7, level=c(80,90))
## h=12 para predecir el año 2010
## Plot sub-meter 3 forecast, limit y and add labels
plot(forecastfitSM2_Semanal, ylab= "Watt-Hours", xlab="Fecha",
     main="Predicción del consumo energético del año 2010 \n  a través de un modelo de regresión")

Se prevee una bajada del consumo

SM3. AC y termo

Modelo

fitSM3_Semanal <- tslm(tsSM3_Semanal ~ trend + season) 
summary(fitSM3_Semanal)

Call:
tslm(formula = tsSM3_Semanal ~ trend + season)

Residuals:
   Min     1Q Median     3Q    Max 
-50070 -20032  -1703  18539  49320 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   410909      20033  20.511 2.76e-11 ***
trend           5485       1171   4.684 0.000427 ***
season2         4110      25070   0.164 0.872309    
season3        -5833      25152  -0.232 0.820231    
season4       -45960      25288  -1.817 0.092260 .  
season5        16778      25477   0.659 0.521682    
season6        36018      25718   1.401 0.184776    
season7       -48786      26010  -1.876 0.083330 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 30670 on 13 degrees of freedom
Multiple R-squared:  0.7601,    Adjusted R-squared:  0.6309 
F-statistic: 5.884 on 7 and 13 DF,  p-value: 0.00305

Resultados:

  • Valor del coeficiente de determinación: 0.7601. Es un valor que podría ser aceptable, el modelo únicamente explica el 76% de la variabilidad total.

  • Hay tres coeficientes significativamente no nulos.

Predicción IC

Predicción para los próximos 12 meses (año 2010)

forecastfitSM3_Semanal <- forecast(fitSM3_Semanal, h=7, level=c(80,90))
## h=12 para predecir el año 2010
## Plot sub-meter 3 forecast, limit y and add labels
plot(forecastfitSM3_Semanal, ylab= "Watt-Hours", xlab="Fecha",
     main="Predicción del consumo energético del año 2010 \n  a través de un modelo de regresión")

Se prevee una bajada del consumo

Comparación de los coeficientes y resultados de cada modelo

DatosRealesSeriesSemanales <- Granularidad_semanas %>% filter(year == 2010 )

RMSE_SM1_GranSemanal<-accuracy(forecastfitSM1_Semanal  ,DatosRealesSeriesSemanales$Sub_metering_1)[4]
  R2_SM1_GranSemanal<-0.975
MASE_SM1_GranSemanal<-accuracy(forecastfitSM1_Semanal  ,DatosRealesSeriesSemanales$Sub_metering_1)[12]
RMSE_SM2_GranSemanal<-accuracy(forecastfitSM2_Semanal  ,DatosRealesSeriesSemanales$Sub_metering_2)[4]
  R2_SM2_GranSemanal<-0.7717
MASE_SM2_GranSemanal<-accuracy(forecastfitSM2_Semanal  ,DatosRealesSeriesSemanales$Sub_metering_2)[12]
RMSE_SM3_GranSemanal<-accuracy(forecastfitSM3_Semanal  ,DatosRealesSeriesSemanales$Sub_metering_3)[4]
  R2_SM3_GranSemanal<-0.7601
MASE_SM3_GranSemanal<-accuracy(forecastfitSM3_Semanal  ,DatosRealesSeriesSemanales$Sub_metering_3)[12]

ResEvolSemanal<-
cbind.data.frame(
      Submedidor=c("Cocina","Lavadero","AC y TermoE"),
      RMSE = c(RMSE_SM1_GranSemanal,
               RMSE_SM2_GranSemanal,
               RMSE_SM3_GranSemanal),
     # MASE = c(MASE_SM1_GranMensual,MASE_SM2_GranMensual,MASE_SM3_GranMensual),
      R2 = c(R2_SM1_GranSemanal,
             R2_SM2_GranSemanal,
             R2_SM3_GranSemanal)
)

ResEvolSemanal %>% kable(booktabs=TRUE) %>% kable_styling(latex_options = "striped")
Submedidor RMSE R2
Cocina 24509.62 0.9750
Lavadero 22414.30 0.7717
AC y TermoE 65773.10 0.7601
# MASE: Mean absolute error
# RMSE: root mean scuare error. Diferencia entre los valores predichos y los observados.

El RMSE del modelo con mejor coeficiente de determinación es el más alto. Pero no hay que olvidar que el consumo eléctrico para el submedidor correspondiente al AC y termo eléctrico es el que más energía consume de todos con una diferencia muy grande.

Forecasting descomponiendo la serie (así ok)

Descomponer: quitar tendencia

Descomposición de la serie y visualización

Submetering 1

 tsSM1_Semanal %>%
  stl( #t.window=13,
      s.window="periodic", 
      robust=TRUE)  %>% summary()
 Call:
 stl(x = ., s.window = "periodic", robust = TRUE)

 Time.series components:
    seasonal             trend            remainder         
 Min.   :-27570.39   Min.   :85332.42   Min.   :-17842.000  
 1st Qu.:-27284.26   1st Qu.:86567.92   1st Qu.: -2325.906  
 Median :-17054.28   Median :88045.73   Median :  -446.063  
 Mean   :     0.00   Mean   :88659.80   Mean   : -1989.231  
 3rd Qu.: 45539.85   3rd Qu.:90927.75   3rd Qu.:   621.535  
 Max.   : 59693.44   Max.   :91708.02   Max.   :  5531.746  
 IQR:
     STL.seasonal STL.trend STL.remainder data 
     72824         4360      2947         65904
   % 110.5          6.6       4.5         100.0

 Weights:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.0000  0.7882  0.9452  0.7781  0.9947  0.9999 

 Other components: List of 5
 $ win  : Named num [1:3] 211 11 7
 $ deg  : Named int [1:3] 0 1 1
 $ jump : Named num [1:3] 22 2 1
 $ inner: int 1
 $ outer: int 15
tsSM1_Semanal %>%
  stl(t.window=13, s.window="periodic", robust=TRUE) %>%
  autoplot()

Remainder no es un muelle. Observamos tendencia creciente y fuerte componente estacional

Se aprecia componente estacional: oscilacione que se repiten cada año. La componente tendencia tambien se aprecia. Tendencia decreciente.

Componentes_SM1_SemanalSTL<- tsSM1_Semanal %>%stl( #t.window=13,
  s.window="periodic",  robust=TRUE)
autoplot(Componentes_SM1_SemanalSTL, ts.colour = 'blue')

plot.ts(tsSM1_Semanal, col = 'gray')
lines(Componentes_SM1_SemanalSTL$time.series[,2], 
      col = "red", lwd = 1, lty = 2)

Submetering 2

 tsSM2_Semanal %>%
  stl( #t.window=13,
      s.window="periodic", 
      robust=TRUE) %>% summary()
 Call:
 stl(x = ., s.window = "periodic", robust = TRUE)

 Time.series components:
    seasonal             trend             remainder        
 Min.   :-39520.48   Min.   : 76016.18   Min.   :-20665.91  
 1st Qu.:-36421.48   1st Qu.: 78225.80   1st Qu.: -5262.30  
 Median : -4900.90   Median : 85975.97   Median :  -975.05  
 Mean   :     0.00   Mean   : 95499.46   Mean   :  5065.54  
 3rd Qu.: 29057.65   3rd Qu.:112604.73   3rd Qu.:  7002.00  
 Max.   : 47277.29   Max.   :133782.17   Max.   : 62291.18  
 IQR:
     STL.seasonal STL.trend STL.remainder data 
     65479        34379     12264         60848
   % 107.6         56.5      20.2         100.0

 Weights:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.0000  0.8936  0.9452  0.8378  0.9882  0.9991 

 Other components: List of 5
 $ win  : Named num [1:3] 211 11 7
 $ deg  : Named int [1:3] 0 1 1
 $ jump : Named num [1:3] 22 2 1
 $ inner: int 1
 $ outer: int 15
tsSM2_Semanal %>%
  stl(t.window=13, s.window="periodic", robust=TRUE) %>%
  autoplot()

Componentes_SM2_SemanalSTL<- tsSM2_Semanal %>%stl( #t.window=13,
  s.window="periodic",  robust=TRUE)
autoplot(Componentes_SM2_SemanalSTL, ts.colour = 'blue')

plot.ts(tsSM2_Semanal, col = 'gray')
lines(Componentes_SM2_SemanalSTL$time.series[,2], 
      col = "red", lwd = 1, lty = 2)

Se aprecia componente estacional: oscilacione que se repiten cada año. La componente tendencia tambien se aprecia. Tendencia decreciente.

Submetering 3

Se aprecia componente estacional: oscilacione que se repiten cada año. La componente tendencia tambien se aprecia. Tendencia creciente

tsSM3_Semanal %>%
  stl( #t.window=13,
      s.window="periodic", 
      robust=TRUE)  %>% summary()
 Call:
 stl(x = ., s.window = "periodic", robust = TRUE)

 Time.series components:
    seasonal             trend            remainder         
 Min.   :-46280.96   Min.   :444191.6   Min.   :-117714.08  
 1st Qu.:-44234.04   1st Qu.:449466.2   1st Qu.: -10712.44  
 Median : -5540.13   Median :456724.5   Median :  -6746.47  
 Mean   :     0.00   Mean   :470965.9   Mean   :  -5957.36  
 3rd Qu.: 40335.65   3rd Qu.:489708.8   3rd Qu.:  11362.30  
 Max.   : 41697.85   Max.   :534471.1   Max.   :  29107.15  
 IQR:
     STL.seasonal STL.trend STL.remainder data 
     84570        40243     22075         63070
   % 134.1         63.8      35.0         100.0

 Weights:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.0000  0.9035  0.9452  0.8755  0.9737  0.9937 

 Other components: List of 5
 $ win  : Named num [1:3] 211 11 7
 $ deg  : Named int [1:3] 0 1 1
 $ jump : Named num [1:3] 22 2 1
 $ inner: int 1
 $ outer: int 15
Componentes_SM3_SemanalSTL<- tsSM3_Semanal %>%stl( #t.window=13,
  s.window="periodic",  robust=TRUE)
autoplot(Componentes_SM3_SemanalSTL, ts.colour = 'blue')

plot.ts(tsSM3_Semanal, col = 'gray')
lines(Componentes_SM3_SemanalSTL$time.series[,2], 
      col = "red", lwd = 1, lty = 2)

Tendencia creciente del consumo energético

Predicción con Holt-Winters (suavizado exponencial) para el año 2010

Sin estacionalidad

Sería interesante verlo también con estacionalidad, para ver si se mantiene o no la tendencia.

Tendencia de predicción (crecimiento o decrecimiento del consumo) y predicción teniendo en cuenta las componentes

Evolución mensual del consumo (sin estacionalidad). Años 2007-2009

Cocina
#tsSM1_Mensual_Adjusted <- tsSM1_Mensual - Componentes_SM1_Mensual$seasonal
#autoplot(tsSM1_Mensual_Adjusted)
tsSM1_Semanal_detrended <- tsSM1_Semanal - Componentes_SM1_SemanalSTL$time.series[,2]
plot.ts(tsSM1_Semanal_detrended, 
        ylab = "Energía (Vatios-hora)", 
        main = 'Seasonal variability of energy comsumption',
        cex.main = 0.85)

par(mfrow = c(1,2))
plot(Componentes_SM1_SemanalSTL$time.series[,3], 
     col = "blue", main = 'Remainder', ylab = "")
qqnorm(Componentes_SM1_SemanalSTL$time.series[,3])
qqline(Componentes_SM1_SemanalSTL$time.series[,3])

shapiro.test(Componentes_SM1_SemanalSTL$time.series[,3])

    Shapiro-Wilk normality test

data:  Componentes_SM1_SemanalSTL$time.series[, 3]
W = 0.80542, p-value = 0.0007949

Los resíduos no son solo ruido. Aún hay estacionalidad

Volvemos a descomponer la serie desestacionalizada una vez.

## Volvemos a descomponer la serie
plot(stl(tsSM1_Semanal_detrended ,s.window = "periodic",robust = TRUE))

La componente estacional se mueve en un ranfo de \(10^{-11}\), podemos asumir que se ha eliminado esta componente.

Suavizado exponencial de HoltWinters
## Holt Winters Exponential Smoothing & Plot
tsSM1_HW_Semanal <- HoltWinters(tsSM1_Semanal_detrended ,
                              #  tsSM1_Mensual_Adjusted, 
  beta=TRUE, gamma=TRUE)
# optimiza
plot(tsSM1_HW_Semanal,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Pronóstico de HoltWinters para el año 2010

Vamos a predecir las próximas 12 observaciones.

## HoltWinters forecast & plot
tsSM1_HW_Semanal_for <- forecast(tsSM1_HW_Semanal, h=7)
plot(tsSM1_HW_Semanal_for , ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico",
     main="Predicción del consumo energético semanal de la cocina \n Año 2010" ) 

## Forecast HoltWinters with diminished confidence levels
tsSM1_HW_Semanal_forC <- forecast(tsSM1_HW_Semanal, h=7, level=c(10,25))
## Plot only the forecasted area
plot(tsSM1_HW_Semanal_forC, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico", main="Predicción del consumo energético semanal para el año 2010", start(2010))

Lavadero
## Seasonal adjusting sub-meter 3 by subtracting the seasonal component & plot
# tsSM2_Mensual_Adjusted <- tsSM2_Mensual - Componentes_SM2_Mensual$seasonal
# autoplot(tsSM2_Mensual_Adjusted)
tsSM2_Semanal_detrended <- tsSM2_Semanal - Componentes_SM2_SemanalSTL$time.series[,2]
plot.ts(tsSM2_Semanal_detrended, 
        ylab = "Energía (Vatios-hora)", 
        main = 'Seasonal variability of energy comsumption',
        cex.main = 0.85)

par(mfrow = c(1,2))
plot(Componentes_SM2_SemanalSTL$time.series[,3], 
     col = "blue", main = 'Remainder', ylab = "")
qqnorm(Componentes_SM2_SemanalSTL$time.series[,3])
qqline(Componentes_SM2_SemanalSTL$time.series[,3])

shapiro.test(Componentes_SM2_SemanalSTL$time.series[,3])

    Shapiro-Wilk normality test

data:  Componentes_SM2_SemanalSTL$time.series[, 3]
W = 0.70675, p-value = 3.345e-05

Los resíduos no siguen una distribución normal. La descomposición no es buena.

## Volvemos a descomponer la serie
plot(stl(tsSM2_Semanal_detrended ,s.window = "periodic",robust = TRUE))

La componente estacional se mueve en un rango de \(10^{-11}\), podemos asumir que se ha eliminado esta componente.

Suavizado exponencial de HoltWinters
## Holt Winters Exponential Smoothing & Plot
tsSM2_HW_Semanal <- HoltWinters(tsSM2_Semanal_detrended ,
                                # tsSM2_Mensual_Adjusted,
                                beta=TRUE, gamma=TRUE)
plot(tsSM2_HW_Semanal,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Pronóstico de HoltWinters para el año 2010

Vamos a predecir las próximas 7 observaciones.

## HoltWinters forecast & plot
     tsSM2_HW_Semanal_for <- forecast(tsSM2_HW_Semanal, h=7)
plot(tsSM2_HW_Semanal_for , ylab= "Vatios-Hora", xlab="Fecha - Lavadero",
     main="Predicción semanal del consumo energético del lavadero \ Año 2010" ) 

## Forecast HoltWinters with diminished confidence levels
tsSM2_HW_Semanal_forC <- forecast(tsSM2_HW_Semanal, h=12, level=c(10,25))
## Plot only the forecasted area
plot(tsSM2_HW_Semanal_forC, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico", main="Predicción semanal para el año 2010", start(2010))

AC y termo eléctrico
# ## Seasonal adjusting sub-meter 3 by subtracting the seasonal component & plot
# tsSM3_Mensual_Adjusted <- tsSM3_Mensual - Componentes_SM3_Mensual$seasonal
# autoplot(tsSM3_Mensual_Adjusted)
tsSM3_Semanal_detrended <- tsSM3_Semanal - Componentes_SM3_SemanalSTL$time.series[,2]
plot.ts(tsSM3_Semanal_detrended, 
        ylab = "Energía (Vatios-hora)", 
        main = 'Seasonal variability of energy comsumption',
        cex.main = 0.85)

par(mfrow = c(1,2))
plot(Componentes_SM3_SemanalSTL$time.series[,3], 
     col = "blue", main = 'Remainder', ylab = "")
qqnorm(Componentes_SM3_SemanalSTL$time.series[,3])
qqline(Componentes_SM3_SemanalSTL$time.series[,3])

shapiro.test(Componentes_SM3_SemanalSTL$time.series[,3])

    Shapiro-Wilk normality test

data:  Componentes_SM3_SemanalSTL$time.series[, 3]
W = 0.67378, p-value = 1.32e-05

Los resíduos no siguen una distribución normal. La descomposición no es buena.

## Volvemos a descomponer la serie
plot(stl(tsSM3_Semanal_detrended ,s.window = "periodic",robust = TRUE))

La componente estacional se mueve en un ranfo de \(10^{-11}\), podemos asumir que se ha eliminado esta componente.

Suavizado exponencial de HoltWinters
## Holt Winters Exponential Smoothing & Plot
tsSM3_HW_Semanal <- HoltWinters(tsSM3_Semanal_detrended ,
                              # tsSM3_Mensual_Adjusted,
                                beta=TRUE, gamma=TRUE)
plot(tsSM3_HW_Semanal,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Pronóstico de HoltWinters para el año 2010

Vamos a predecir las próximas 12 observaciones.

## HoltWinters forecast & plot
tsSM3_HW_Semanal_for <- forecast(tsSM3_HW_Semanal, h=7)
plot(tsSM3_HW_Semanal_for , ylab= "Vatios-Hora", xlab="Fecha - Lavadero",
     main="Predicción del consumo energético semanal del AC y termo eléctrico \n Año 2010" ) 

## Forecast HoltWinters with diminished confidence levels
tsSM3_HW_Semanal_forC <- forecast(tsSM3_HW_Semanal, h=7, level=c(10,25))
## Plot only the forecasted area
plot(tsSM3_HW_Semanal_forC, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico", main="Predicción semanal del consumo energético para el año 2010", start(2010))

Comparación de los coeficientes y resultados de cada predicción. Datos sin estacionalidad

En las gráficas de predicción, hemos observado que el ajuste no mantiene la tendencia de consumo energético de la serie.

RMSE_SM1_GranSemanalFor<-accuracy(tsSM1_HW_Semanal_forC   ,DatosRealesSeriesSemanales$Sub_metering_1)[4]
RMSE_SM2_GranSemanalFor<-accuracy(tsSM2_HW_Semanal_forC   ,DatosRealesSeriesSemanales$Sub_metering_2)[4]
RMSE_SM3_GranSemanalFor<-accuracy(tsSM3_HW_Semanal_forC   ,DatosRealesSeriesSemanales$Sub_metering_3)[4]


ResEvolSemanal2_SinEstacionalidadConTendencia <-
cbind.data.frame(
      Submedidor=c("Cocina","Lavadero","AC y TermoE"),
      RMSE = c(RMSE_SM1_GranSemanalFor,
               RMSE_SM2_GranSemanalFor,
               RMSE_SM3_GranSemanalFor)
)
ResEvolSemanal2_SinEstacionalidadConTendencia %>% kable(booktabs=TRUE) %>% kable_styling(latex_options = "striped")
Submedidor RMSE
Cocina 59793.12
Lavadero 63672.21
AC y TermoE 481966.03

Evolución mensual del consumo (con estacionalidad). Años 2007-2009

Vamos a ver una predicción del consumo enegético de cada submedidor con el método de HW. En este caso, no eliminaremos la estacionalidad, para ver así si se mantiene la tendencia de consumo.

Cocina
tsSM1_Mensual_AdjustedII <- tsSM1_Mensual # - Componentes_SM1_Mensual$seasonal
autoplot(tsSM1_Mensual_AdjustedII)

Parece un muelle.

## Volvemos a descomponer la serie
plot(decompose(tsSM1_Mensual_AdjustedII))

La componente estacional se mueve en un rango de \(10^{-11}\), podemos asumir que se ha eliminado esta componente.

Suavizado exponencial de HoltWinters
## Holt Winters Exponential Smoothing & Plot
tsSM1_HW_MensualII <- HoltWinters(tsSM1_Mensual_AdjustedII, beta=TRUE, gamma=TRUE)
# optimiza
plot(tsSM1_HW_MensualII,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Pronóstico de HoltWinters para el año 2010

Vamos a predecir las próximas 12 observaciones.

## HoltWinters forecast & plot
tsSM1_HW_Mensual_forII <- forecast(tsSM1_HW_MensualII, h=12)
plot(tsSM1_HW_Mensual_forII , ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico",
     main="Predicción del consumo energético de la cocina \ Año 2010" ) 

## Forecast HoltWinters with diminished confidence levels
tsSM1_HW_Mensual_forCII <- forecast(tsSM1_HW_MensualII, h=12, level=c(10,25))
## Plot only the forecasted area
plot(tsSM1_HW_Mensual_forCII, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico", main="Predicción para el año 2010", start(2010))

Lavadero
## Seasonal adjusting sub-meter 3 by subtracting the seasonal component & plot
tsSM2_Mensual_AdjustedII <- tsSM2_Mensual ## - Componentes_SM2_Mensual$seasonal
autoplot(tsSM2_Mensual_AdjustedII)

No arece un muelle.

## Volvemos a descomponer la serie
plot(decompose(tsSM2_Mensual_AdjustedII)) 

TENDENCIA CLARAMENTE DECRECIENTE. Tambien observamos esacionalidad

Suavizado exponencial de HoltWinters
## Holt Winters Exponential Smoothing & Plot
tsSM2_HW_MensualII <- HoltWinters(tsSM2_Mensual_AdjustedII, beta=TRUE, gamma=TRUE)
plot(tsSM2_HW_MensualII,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Pronóstico de HoltWinters para el año 2010

Vamos a predecir las próximas 12 observaciones.

## HoltWinters forecast & plot
tsSM2_HW_Mensual_forII <- forecast(tsSM2_HW_MensualII, h=12)
plot(tsSM2_HW_Mensual_forII , ylab= "Vatios-Hora", xlab="Fecha - Lavadero",
     main="Predicción del consumo energético del lavadero \ Año 2010" ) 

## Forecast HoltWinters with diminished confidence levels
tsSM2_HW_Mensual_forCII <- forecast(tsSM2_HW_MensualII, h=12, level=c(10,25))
## Plot only the forecasted area
plot(tsSM2_HW_Mensual_forCII, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico", main="Predicción para el año 2010", start(2010))

AC y termo eléctrico
## Seasonal adjusting sub-meter 3 by subtracting the seasonal component & plot
tsSM3_Mensual_AdjustedII <- tsSM3_Mensual # - Componentes_SM3_Mensual$seasonal
autoplot(tsSM3_Mensual_AdjustedII)

No parece un muelle.

## Volvemos a descomponer la serie
plot(decompose(tsSM3_Mensual_AdjustedII))

Suavizado exponencial de HoltWinters
## Holt Winters Exponential Smoothing & Plot
tsSM3_HW_MensualII <- HoltWinters(tsSM3_Mensual_AdjustedII, beta=TRUE, gamma=TRUE)
plot(tsSM3_HW_MensualII,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Pronóstico de HoltWinters para el año 2010

Vamos a predecir las próximas 12 observaciones.

## HoltWinters forecast & plot
tsSM3_HW_Mensual_forII <- forecast(tsSM3_HW_MensualII, h=12)
plot(tsSM3_HW_Mensual_forII , ylab= "Vatios-Hora", xlab="Fecha - Lavadero",
     main="Predicción del consumo energético del AC y termo eléctrico \ Año 2010" ) 

## Forecast HoltWinters with diminished confidence levels
tsSM3_HW_Mensual_forCII <- forecast(tsSM3_HW_MensualII, h=12, level=c(10,25))
## Plot only the forecasted area
plot(tsSM3_HW_Mensual_forCII, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico", main="Predicción para el año 2010", start(2010))

Comparación de los coeficientes y resultados de cada predicción. Datos con estacionalidad
RMSE_SM1_GranMensualForII<-accuracy(tsSM1_HW_Mensual_forCII   ,DatosRealesSeriesMensuales$Sub_metering_1)[4]
RMSE_SM2_GranMensualForII<-accuracy(tsSM2_HW_Mensual_forCII   ,DatosRealesSeriesMensuales$Sub_metering_2)[4]
RMSE_SM3_GranMensualForII<-accuracy(tsSM3_HW_Mensual_forCII   ,DatosRealesSeriesMensuales$Sub_metering_3)[4]


ResEvolMensual2SinTendenciaConEstacionalidad<-
cbind.data.frame(
      Submedidor=c("Cocina","Lavadero","AC y TermoE"),
      RMSE = c(RMSE_SM1_GranMensualForII,
               RMSE_SM2_GranMensualForII,
               RMSE_SM3_GranMensualForII)
   
    
)

ResEvolMensual2SinTendenciaConEstacionalidad %>% kable(booktabs=TRUE) %>% kable_styling(latex_options = "striped")
Submedidor RMSE
Cocina 16362.16
Lavadero 27223.90
AC y TermoE 165249.73

Comparación de resultados. Predicción antes y despues de eliminar estacionalidad

Comparacion<-cbind.data.frame(
  ResEvolMensual2_SinEstacionalidadConTendencia,
  ResEvolMensual2SinTendenciaConEstacionalidad[,2])

colnames(Comparacion) <- c("Submedidor","RMSE Con Estacionalidad","RMSE Sin estacionalidad")
Comparacion %>% kable(booktabs=TRUE) %>% kable_styling(latex_options = "striped")
Submedidor RMSE Con Estacionalidad RMSE Sin estacionalidad
Cocina 45122.42 16362.16
Lavadero 44294.27 27223.90
AC y TermoE 318158.49 165249.73

Resultados más fiables y con menor error trás haber eliminado la estacionalidad de la serie. También así podemos ver la tendencia.

#p1<-autoplot(tsSM1_Mensual_Adjusted)
#p2<-autoplot(tsSM1_Mensual_AdjustedII)
#grid.arrange(
#             grobs = list(p1,p2),
#           nrow = 2,
#             top="Sin estacionalidad",bottom="Con estacionalidad"
#             )

Evolución del consumo energético por semana del año.

Serie temporal y visualización

Vamos a representar la serie por meses para cada submedidor y tambien para la energía global co granularidad mensual. En este caso, la frecuencia será 12, tenemos una observación por mes para los años 2007,8 y 2009.

Granularidad_semanas <- Granularidad_horas %>% group_by(year, week) %>% 
   summarize(
      Sub_metering_1 = sum (Sub_metering_1),
      Sub_metering_2 = sum (Sub_metering_2),
      Sub_metering_3 = sum (Sub_metering_3),
      Global_active_power = sum (Global_active_power),
      energia2 = sum (energia2)
   )
seriesSemanales <- rbind.data.frame(
   filter(Granularidad_semanas, year== 2007  ),
   filter(Granularidad_semanas, year== 2008  ),
   filter(Granularidad_semanas, year== 2009  )
)

## Creamos un objeto TS
   tsSM1_Semanal<- ts(seriesSemanales$Sub_metering_1,         frequency=53, start=c(2007,1))
   tsSM2_Semanal<- ts(seriesSemanales$Sub_metering_2,         frequency=53, start=c(2007,1))
   tsSM3_Semanal<- ts(seriesSemanales$Sub_metering_3,         frequency=53, start=c(2007,1))
tsSM_GAP_Semanal<- ts(seriesSemanales$Global_active_power,    frequency=53, start=c(2007,1))
## Plot sub-meter 3 with autoplot - add labels, color
autoplot(tsSM1_Semanal, ts.colour = 'red', xlab = "Semana del año", ylab = "Energía (Varios-hora)", main = "Evolución semanal del consumo de energía de la cocina \n Años 2007-2009")

autoplot(tsSM2_Semanal, xlab = "Semana del año", ylab = "Energía (Varios-hora)", main = "Evolución semanal del consumo de energía de la lavandería \n Años 2007-2009")

autoplot(tsSM3_Semanal, ts.colour = 'blue', xlab = "Semana del año", ylab = "Energía (Varios-hora)", main = "Evolución semanal del consumo de energía del aire acondicionado y termo eléctrico \n Años 2007-2009")

autoplot(tsSM_GAP_Semanal, ts.colour = 'green', xlab = "Semana del año", ylab = "Energía (Varios-hora)", main = "Evolución semanal del consumo de energía \n Años 2007-2009")

AC y Termo: podría apreciarse cierta estacionalidad Energía total: muelle

Predicción del año 2010 antes de descomponer la serie

Vamos a plicar regresión lineal a cada serie temporal, y veremos el valor de \(R^2\) y del error cuadrático medio.

SM1. Cocina

Modelo

fitSM1_Semanal <- tslm(tsSM1_Semanal ~ trend + season) 
summary(fitSM1_Semanal)

Call:
tslm(formula = tsSM1_Semanal ~ trend + season)

Residuals:
     Min       1Q   Median       3Q      Max 
-10942.7  -1970.9   -182.3   2205.3  12203.0 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 10835.442   2693.122   4.023 0.000108 ***
trend          -8.872      8.426  -1.053 0.294745    
season2      6872.206   3753.914   1.831 0.069985 .  
season3      6643.745   3753.942   1.770 0.079662 .  
season4      4273.617   3753.989   1.138 0.257536    
season5      3290.489   3754.056   0.877 0.382751    
season6      4326.029   3754.141   1.152 0.251800    
season7       791.901   3754.245   0.211 0.833347    
season8      -313.227   3754.368  -0.083 0.933669    
season9     -3984.688   3754.509  -1.061 0.290985    
season10     6229.851   3754.670   1.659 0.100054    
season11     2045.390   3754.850   0.545 0.587092    
season12     5559.929   3755.048   1.481 0.141693    
season13     3915.802   3755.266   1.043 0.299459    
season14     1586.341   3755.502   0.422 0.673594    
season15     3100.547   3755.757   0.826 0.410934    
 [ reached getOption("max.print") -- omitted 38 rows ]
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4598 on 105 degrees of freedom
Multiple R-squared:  0.4836,    Adjusted R-squared:  0.223 
F-statistic: 1.855 on 53 and 105 DF,  p-value: 0.003663

Resultados:

  • Valor del coeficiente de determinación: 0.4836. Es un valor bajo, el modelo únicamente explica el 48% de la variabilidad total.

  • Hay un coeficiente significativamente no nulo.

Predicción IC

Predicción para los próximos 20 días

forecastfitSM1_Semanal <- forecast(fitSM1_Semanal, h=53, level=c(80,90))
## h=12 para predecir el año 2010
## Plot sub-meter 3 forecast, limit y and add labels
plot(forecastfitSM1_Semanal, ylab= "Watt-Hours", xlab="Fecha",
     main="Predicción del consumo energético semanal total del año 2010 \n  a través de un modelo de regresión")

SM2. Lavadero

Modelo

fitSM2_Semanal <- tslm(tsSM2_Semanal ~ trend + season) 
summary(fitSM2_Semanal)

Call:
tslm(formula = tsSM2_Semanal ~ trend + season)

Residuals:
   Min     1Q Median     3Q    Max 
-12710  -2492    262   2070   9448 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  19974.287   2851.564   7.005 2.45e-10 ***
trend          -46.079      8.921  -5.165 1.15e-06 ***
season2      -4864.254   3974.765  -1.224 0.223773    
season3        714.159   3974.795   0.180 0.857756    
season4       -986.429   3974.845  -0.248 0.804489    
season5      -1733.682   3974.915  -0.436 0.663619    
season6       -145.270   3975.005  -0.037 0.970917    
season7      -4670.524   3975.115  -1.175 0.242676    
season8      -1727.111   3975.245  -0.434 0.664841    
season9      -5223.698   3975.395  -1.314 0.191707    
season10      3920.048   3975.565   0.986 0.326383    
season11      1273.127   3975.756   0.320 0.749436    
season12     -1349.793   3975.966  -0.339 0.734920    
season13      -436.381   3976.196  -0.110 0.912818    
season14     -2924.635   3976.446  -0.735 0.463681    
season15     -2526.222   3976.716  -0.635 0.526645    
 [ reached getOption("max.print") -- omitted 38 rows ]
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4868 on 105 degrees of freedom
Multiple R-squared:  0.5326,    Adjusted R-squared:  0.2966 
F-statistic: 2.257 on 53 and 105 DF,  p-value: 0.0001984

Resultados:

  • Valor del coeficiente de determinación: 0.5326. Es un valor aceptable, el modelo únicamente explica el 53% de la variabilidad total.

  • Hay tres coeficientes significativamente no nulos.

Predicción IC

Predicción para los próximos 12 meses (año 2010)

forecastfitSM2_Semanal <- forecast(fitSM2_Semanal, h=53, level=c(80,90))
## h=12 para predecir el año 2010
## Plot sub-meter 3 forecast, limit y and add labels
plot(forecastfitSM2_Semanal, ylab= "Watt-Hours", xlab="Fecha",
     main="Predicción del consumo energético del año 2010 \n  a través de un modelo de regresión")

Se prevee una bajada del consumo

SM3. AC y termo

Modelo

fitSM3_Semanal <- tslm(tsSM3_Semanal ~ trend + season) 
summary(fitSM3_Semanal)

Call:
tslm(formula = tsSM3_Semanal ~ trend + season)

Residuals:
   Min     1Q Median     3Q    Max 
-35933  -5644    -17   6305  26151 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  56824.99    7005.03   8.112 9.94e-13 ***
trend           95.69      21.92   4.366 2.97e-05 ***
season2      13899.65    9764.23   1.424  0.15755    
season3      18434.30    9764.31   1.888  0.06180 .  
season4       9370.94    9764.43   0.960  0.33941    
season5      20968.59    9764.60   2.147  0.03406 *  
season6      14556.57    9764.82   1.491  0.13903    
season7       8915.22    9765.09   0.913  0.36335    
season8      -2728.46    9765.41  -0.279  0.78049    
season9     -25352.82    9765.78  -2.596  0.01078 *  
season10      2945.17    9766.20   0.302  0.76358    
season11      8872.81    9766.67   0.908  0.36571    
season12      8371.46    9767.18   0.857  0.39334    
season13     13322.78    9767.75   1.364  0.17550    
season14      9586.09    9768.36   0.981  0.32868    
season15     -1125.59    9769.03  -0.115  0.90849    
 [ reached getOption("max.print") -- omitted 38 rows ]
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 11960 on 105 degrees of freedom
Multiple R-squared:  0.7153,    Adjusted R-squared:  0.5715 
F-statistic: 4.977 on 53 and 105 DF,  p-value: 1.341e-12

Resultados:

  • Valor del coeficiente de determinación: 0.7153. Es un valor que podría ser aceptable, el modelo únicamente explica el 71% de la variabilidad total.

  • Hay tres coeficientes significativamente no nulos.

Predicción IC

Predicción para los próximos 12 meses (año 2010)

forecastfitSM3_Semanal <- forecast(fitSM3_Semanal, h=53, level=c(80,90))
## h=12 para predecir el año 2010
## Plot sub-meter 3 forecast, limit y and add labels
plot(forecastfitSM3_Semanal, ylab= "Watt-Hours", xlab="Fecha",
     main="Predicción del consumo energético del año 2010 \n  a través de un modelo de regresión")

Se prevee una bajada del consumo

Comparación de los coeficientes y resultados de cada modelo

###  DatosRealesSeriesSemanales <- Granularidad_semanas %>% filter(year == 2010 ###  )
###  
###  RMSE_SM1_GranSemanal<-accuracy(forecastfitSM1_Semanal  ###  ,DatosRealesSeriesSemanales$Sub_metering_1)[4]
###    R2_SM1_GranSemanal<-0.975
###  MASE_SM1_GranSemanal<-accuracy(forecastfitSM1_Semanal  ###  ,DatosRealesSeriesSemanales$Sub_metering_1)[12]
###  RMSE_SM2_GranSemanal<-accuracy(forecastfitSM2_Semanal  ###  ,DatosRealesSeriesSemanales$Sub_metering_2)[4]
###    R2_SM2_GranSemanal<-0.7717
###  MASE_SM2_GranSemanal<-accuracy(forecastfitSM2_Semanal  ###  ,DatosRealesSeriesSemanales$Sub_metering_2)[12]
###  RMSE_SM3_GranSemanal<-accuracy(forecastfitSM3_Semanal  ###  ,DatosRealesSeriesSemanales$Sub_metering_3)[4]
###    R2_SM3_GranSemanal<-0.7601
###  MASE_SM3_GranSemanal<-accuracy(forecastfitSM3_Semanal  ###  ,DatosRealesSeriesSemanales$Sub_metering_3)[12]
###  
###  ResEvolSemanal<-
###  cbind.data.frame(
###        Submedidor=c("Cocina","Lavadero","AC y TermoE"),
###        RMSE = c(RMSE_SM1_GranSemanal,
###                 RMSE_SM2_GranSemanal,
###                 RMSE_SM3_GranSemanal),
###       # MASE = c(MASE_SM1_GranMensual,MASE_SM2_GranMensual,MASE_SM3_GranMen###  sual),
###        R2 = c(R2_SM1_GranSemanal,
###               R2_SM2_GranSemanal,
###               R2_SM3_GranSemanal))
###  ResEvolSemanal %>% kable(booktabs=TRUE) %>% kable_styling(latex_options = ###  "striped")
###  
###  # MASE: Mean absolute error
###  # RMSE: root mean scuare error. Diferencia entre los valores predichos y ###  los observados.

El RMSE del modelo con mejor coeficiente de determinación es el más alto. Pero no hay que olvidar que el consumo eléctrico para el submedidor correspondiente al AC y termo eléctrico es el que más energía consume de todos con una diferencia muy grande.

Forecasting descomponiendo la serie (así ok)

Descomponer: quitar tendencia

Descomposición de la serie y visualización

Submetering 1

 tsSM1_Semanal %>%
  stl( #t.window=13,
      s.window="periodic", 
      robust=TRUE)  %>% summary()
 Call:
 stl(x = ., s.window = "periodic", robust = TRUE)

 Time.series components:
    seasonal              trend             remainder         
 Min.   :-11380.560   Min.   : 9373.389   Min.   :-17306.164  
 1st Qu.: -1525.527   1st Qu.:10411.060   1st Qu.: -1565.107  
 Median :   605.971   Median :11525.838   Median :    50.446  
 Mean   :     0.000   Mean   :11285.767   Mean   :   161.290  
 3rd Qu.:  2468.518   3rd Qu.:12047.600   3rd Qu.:  1583.123  
 Max.   : 12092.727   Max.   :12804.151   Max.   : 18204.879  
 IQR:
     STL.seasonal STL.trend STL.remainder data 
     3994         1637      3148          6386 
   %  62.5         25.6      49.3         100.0

 Weights:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.0000  0.8303  0.9452  0.8244  0.9824  1.0000 

 Other components: List of 5
 $ win  : Named num [1:3] 1591 81 53
 $ deg  : Named int [1:3] 0 1 1
 $ jump : Named num [1:3] 160 9 6
 $ inner: int 1
 $ outer: int 15
#tsSM1_Semanal %>%
#  stl(t.window=13, s.window="periodic", robust=TRUE) %>%
#  autoplot()

Remainder no es un muelle. Observamos tendencia creciente y fuerte componente estacional

Se aprecia componente estacional: oscilacione que se repiten cada año. La componente tendencia tambien se aprecia. Tendencia decreciente.

Componentes_SM1_SemanalSTL<- tsSM1_Semanal %>%stl( #t.window=13,
  s.window="periodic",  robust=TRUE)
autoplot(Componentes_SM1_SemanalSTL, ts.colour = 'blue')

plot.ts(tsSM1_Semanal, col = 'gray')
lines(Componentes_SM1_SemanalSTL$time.series[,2], 
      col = "red", lwd = 1, lty = 2)

Submetering 2

 tsSM2_Semanal %>%
  stl( #t.window=13,
      s.window="periodic", 
      robust=TRUE) %>% summary()
 Call:
 stl(x = ., s.window = "periodic", robust = TRUE)

 Time.series components:
    seasonal              trend            remainder         
 Min.   :-11082.099   Min.   :10972.07   Min.   :-17744.999  
 1st Qu.: -2312.887   1st Qu.:11343.85   1st Qu.: -1509.064  
 Median :   306.932   Median :12813.79   Median :   128.053  
 Mean   :     0.000   Mean   :13181.15   Mean   :   101.019  
 3rd Qu.:  2737.120   3rd Qu.:15219.46   3rd Qu.:  1570.218  
 Max.   :  9565.068   Max.   :15803.51   Max.   : 14301.230  
 IQR:
     STL.seasonal STL.trend STL.remainder data 
     5050         3876      3079          7916 
   %  63.8         49.0      38.9         100.0

 Weights:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.0000  0.7600  0.9452  0.7848  0.9927  1.0000 

 Other components: List of 5
 $ win  : Named num [1:3] 1591 81 53
 $ deg  : Named int [1:3] 0 1 1
 $ jump : Named num [1:3] 160 9 6
 $ inner: int 1
 $ outer: int 15
#tsSM2_Semanal %>%
#  stl(t.window=13, s.window="periodic", robust=TRUE) %>%
#  autoplot()
Componentes_SM2_SemanalSTL<- tsSM2_Semanal %>%stl( #t.window=13,
  s.window="periodic",  robust=TRUE)
autoplot(Componentes_SM2_SemanalSTL, ts.colour = 'blue')

plot.ts(tsSM2_Semanal, col = 'gray')
lines(Componentes_SM2_SemanalSTL$time.series[,2], 
      col = "red", lwd = 1, lty = 2)

Se aprecia componente estacional: oscilacione que se repiten cada año. La componente tendencia tambien se aprecia. Tendencia decreciente.

Submetering 3

Se aprecia componente estacional: oscilacione que se repiten cada año. La componente tendencia tambien se aprecia. Tendencia creciente

tsSM3_Semanal %>%
  stl( #t.window=13,
      s.window="periodic", 
      robust=TRUE)  %>% summary()
 Call:
 stl(x = ., s.window = "periodic", robust = TRUE)

 Time.series components:
    seasonal             trend            remainder        
 Min.   :-47211.03   Min.   :57753.65   Min.   :-56145.14  
 1st Qu.:-10494.71   1st Qu.:58324.34   1st Qu.: -4929.51  
 Median :  4528.45   Median :63141.11   Median :  -321.52  
 Mean   :     0.00   Mean   :62123.83   Mean   :  -707.60  
 3rd Qu.: 11265.72   3rd Qu.:63648.44   3rd Qu.:  3997.01  
 Max.   : 26621.30   Max.   :68148.68   Max.   : 39973.53  
 IQR:
     STL.seasonal STL.trend STL.remainder data 
     21760         5324      8927         23820
   %  91.4         22.4      37.5         100.0

 Weights:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.0000  0.8598  0.9452  0.8486  0.9898  1.0000 

 Other components: List of 5
 $ win  : Named num [1:3] 1591 81 53
 $ deg  : Named int [1:3] 0 1 1
 $ jump : Named num [1:3] 160 9 6
 $ inner: int 1
 $ outer: int 15
Componentes_SM3_SemanalSTL<- tsSM3_Semanal %>%stl( #t.window=13,
  s.window="periodic",  robust=TRUE)
autoplot(Componentes_SM3_SemanalSTL, ts.colour = 'blue')

plot.ts(tsSM3_Semanal, col = 'gray')
lines(Componentes_SM3_SemanalSTL$time.series[,2], 
      col = "red", lwd = 1, lty = 2)

Tendencia creciente del consumo energético

Predicción con Holt-Winters (suavizado exponencial) para el año 2010

Sin estacionalidad

Sería interesante verlo también con estacionalidad, para ver si se mantiene o no la tendencia.

Tendencia de predicción (crecimiento o decrecimiento del consumo) y predicción teniendo en cuenta las componentes

Evolución mensual del consumo (sin estacionalidad). Años 2007-2009

Cocina
#tsSM1_Mensual_Adjusted <- tsSM1_Mensual - Componentes_SM1_Mensual$seasonal
#autoplot(tsSM1_Mensual_Adjusted)
tsSM1_Semanal_detrended <- tsSM1_Semanal - Componentes_SM1_SemanalSTL$time.series[,2]
plot.ts(tsSM1_Semanal_detrended, 
        ylab = "Energía (Vatios-hora)", 
        main = 'Seasonal variability of energy comsumption',
        cex.main = 0.85)

par(mfrow = c(1,2))
plot(Componentes_SM1_SemanalSTL$time.series[,3], 
     col = "blue", main = 'Remainder', ylab = "")
qqnorm(Componentes_SM1_SemanalSTL$time.series[,3])
qqline(Componentes_SM1_SemanalSTL$time.series[,3])

shapiro.test(Componentes_SM1_SemanalSTL$time.series[,3])

    Shapiro-Wilk normality test

data:  Componentes_SM1_SemanalSTL$time.series[, 3]
W = 0.83235, p-value = 3.231e-12

Los resíduos no son solo ruido. Aún hay estacionalidad

Volvemos a descomponer la serie desestacionalizada una vez.

## Volvemos a descomponer la serie
plot(stl(tsSM1_Semanal_detrended ,s.window = "periodic",robust = TRUE))

La componente estacional se mueve en un ranfo de \(10^{-11}\), podemos asumir que se ha eliminado esta componente.

Suavizado exponencial de HoltWinters
## Holt Winters Exponential Smoothing & Plot
tsSM1_HW_Semanal <- HoltWinters(tsSM1_Semanal_detrended ,
                              #  tsSM1_Mensual_Adjusted, 
  beta=TRUE, gamma=TRUE)
# optimiza
plot(tsSM1_HW_Semanal,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Pronóstico de HoltWinters para el año 2010

Vamos a predecir las próximas 53 observaciones, es decir, el consumo energético semanal de cada submedidor para las 53 semanas del año 2010.

## HoltWinters forecast & plot
tsSM1_HW_Semanal_for <- forecast(tsSM1_HW_Semanal, h=53)
plot(tsSM1_HW_Semanal_for , ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico",
     main="Predicción del consumo energético semanal de la cocina \n Año 2010" ) 

## Forecast HoltWinters with diminished confidence levels
tsSM1_HW_Semanal_forC <- forecast(tsSM1_HW_Semanal, h=53, level=c(10,25))
## Plot only the forecasted area
plot(tsSM1_HW_Semanal_forC, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico", main="Predicción del consumo energético semanal para el año 2010", start(2010))

Lavadero
## Seasonal adjusting sub-meter 3 by subtracting the seasonal component & plot
# tsSM2_Mensual_Adjusted <- tsSM2_Mensual - Componentes_SM2_Mensual$seasonal
# autoplot(tsSM2_Mensual_Adjusted)
tsSM2_Semanal_detrended <- tsSM2_Semanal - Componentes_SM2_SemanalSTL$time.series[,2]
plot.ts(tsSM2_Semanal_detrended, 
        ylab = "Energía (Vatios-hora)", 
        main = 'Seasonal variability of energy comsumption',
        cex.main = 0.85)

par(mfrow = c(1,2))
plot(Componentes_SM2_SemanalSTL$time.series[,3], 
     col = "blue", main = 'Remainder', ylab = "")
qqnorm(Componentes_SM2_SemanalSTL$time.series[,3])
qqline(Componentes_SM2_SemanalSTL$time.series[,3])

shapiro.test(Componentes_SM2_SemanalSTL$time.series[,3])

    Shapiro-Wilk normality test

data:  Componentes_SM2_SemanalSTL$time.series[, 3]
W = 0.88071, p-value = 5.464e-10

Los resíduos no siguen una distribución normal. La descomposición no es buena.

## Volvemos a descomponer la serie
plot(stl(tsSM2_Semanal_detrended ,s.window = "periodic",robust = TRUE))

La componente estacional se mueve en un rango de \(10^{-11}\), podemos asumir que se ha eliminado esta componente.

Suavizado exponencial de HoltWinters
## Holt Winters Exponential Smoothing & Plot
tsSM2_HW_Semanal <- HoltWinters(tsSM2_Semanal_detrended ,
                                # tsSM2_Mensual_Adjusted,
                                beta=TRUE, gamma=TRUE)
plot(tsSM2_HW_Semanal,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Pronóstico de HoltWinters para el año 2010

Vamos a predecir las próximas 7 observaciones.

## HoltWinters forecast & plot
     tsSM2_HW_Semanal_for <- forecast(tsSM2_HW_Semanal, h=53)
plot(tsSM2_HW_Semanal_for , ylab= "Vatios-Hora", xlab="Fecha - Lavadero",
     main="Predicción semanal del consumo energético del lavadero \ Año 2010" ) 

## Forecast HoltWinters with diminished confidence levels
tsSM2_HW_Semanal_forC <- forecast(tsSM2_HW_Semanal, h=53, level=c(10,25))
## Plot only the forecasted area
plot(tsSM2_HW_Semanal_forC, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico", main="Predicción semanal para el año 2010", start(2010))

AC y termo eléctrico
# ## Seasonal adjusting sub-meter 3 by subtracting the seasonal component & plot
# tsSM3_Mensual_Adjusted <- tsSM3_Mensual - Componentes_SM3_Mensual$seasonal
# autoplot(tsSM3_Mensual_Adjusted)
tsSM3_Semanal_detrended <- tsSM3_Semanal - Componentes_SM3_SemanalSTL$time.series[,2]
plot.ts(tsSM3_Semanal_detrended, 
        ylab = "Energía (Vatios-hora)", 
        main = 'Seasonal variability of energy comsumption',
        cex.main = 0.85)

par(mfrow = c(1,2))
plot(Componentes_SM3_SemanalSTL$time.series[,3], 
     col = "blue", main = 'Remainder', ylab = "")
qqnorm(Componentes_SM3_SemanalSTL$time.series[,3])
qqline(Componentes_SM3_SemanalSTL$time.series[,3])

shapiro.test(Componentes_SM3_SemanalSTL$time.series[,3])

    Shapiro-Wilk normality test

data:  Componentes_SM3_SemanalSTL$time.series[, 3]
W = 0.83537, p-value = 4.311e-12

Los resíduos no siguen una distribución normal. La descomposición no es buena.

## Volvemos a descomponer la serie
plot(stl(tsSM3_Semanal_detrended ,s.window = "periodic",robust = TRUE))

La componente estacional se mueve en un ranfo de \(10^{-11}\), podemos asumir que se ha eliminado esta componente.

Suavizado exponencial de HoltWinters
## Holt Winters Exponential Smoothing & Plot
tsSM3_HW_Semanal <- HoltWinters(tsSM3_Semanal_detrended ,
                              # tsSM3_Mensual_Adjusted,
                                beta=TRUE, gamma=TRUE)
plot(tsSM3_HW_Semanal,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Pronóstico de HoltWinters para el año 2010

Vamos a predecir las próximas 12 observaciones.

## HoltWinters forecast & plot
tsSM3_HW_Semanal_for <- forecast(tsSM3_HW_Semanal, h=53)
plot(tsSM3_HW_Semanal_for , ylab= "Vatios-Hora", xlab="Fecha - Lavadero",
     main="Predicción del consumo energético semanal del AC y termo eléctrico \n Año 2010" ) 

## Forecast HoltWinters with diminished confidence levels
tsSM3_HW_Semanal_forC <- forecast(tsSM3_HW_Semanal, h=53, level=c(10,25))
## Plot only the forecasted area
plot(tsSM3_HW_Semanal_forC, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico", main="Predicción semanal del consumo energético para el año 2010", start(2010))

Comparación de los coeficientes y resultados de cada predicción. Datos sin estacionalidad

En las gráficas de predicción, hemos observado que el ajuste no mantiene la tendencia de consumo energético de la serie.

RMSE_SM1_GranSemanalFor<-accuracy(tsSM1_HW_Semanal_forC   ,DatosRealesSeriesSemanales$Sub_metering_1)[4]
RMSE_SM2_GranSemanalFor<-accuracy(tsSM2_HW_Semanal_forC   ,DatosRealesSeriesSemanales$Sub_metering_2)[4]
RMSE_SM3_GranSemanalFor<-accuracy(tsSM3_HW_Semanal_forC   ,DatosRealesSeriesSemanales$Sub_metering_3)[4]


ResEvolSemanal2_SinEstacionalidadConTendencia <-
cbind.data.frame(
      Submedidor=c("Cocina","Lavadero","AC y TermoE"),
      RMSE = c(RMSE_SM1_GranSemanalFor,
               RMSE_SM2_GranSemanalFor,
               RMSE_SM3_GranSemanalFor)
)
ResEvolSemanal2_SinEstacionalidadConTendencia %>% kable(booktabs=TRUE) %>% kable_styling(latex_options = "striped")
Submedidor RMSE
Cocina 64689.54
Lavadero 91145.26
AC y TermoE 545792.05

Evolución mensual del consumo (con estacionalidad). Años 2007-2009

Vamos a ver una predicción del consumo enegético de cada submedidor con el método de HW. En este caso, no eliminaremos la estacionalidad, para ver así si se mantiene la tendencia de consumo.

Cocina
tsSM1_Mensual_AdjustedII <- tsSM1_Mensual # - Componentes_SM1_Mensual$seasonal
autoplot(tsSM1_Mensual_AdjustedII)

Parece un muelle.

## Volvemos a descomponer la serie
plot(decompose(tsSM1_Mensual_AdjustedII))

La componente estacional se mueve en un rango de \(10^{-11}\), podemos asumir que se ha eliminado esta componente.

Suavizado exponencial de HoltWinters
## Holt Winters Exponential Smoothing & Plot
tsSM1_HW_MensualII <- HoltWinters(tsSM1_Mensual_AdjustedII, beta=TRUE, gamma=TRUE)
# optimiza
plot(tsSM1_HW_MensualII,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Pronóstico de HoltWinters para el año 2010

Vamos a predecir las próximas 12 observaciones.

## HoltWinters forecast & plot
tsSM1_HW_Mensual_forII <- forecast(tsSM1_HW_MensualII, h=12)
plot(tsSM1_HW_Mensual_forII , ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico",
     main="Predicción del consumo energético de la cocina \ Año 2010" ) 

## Forecast HoltWinters with diminished confidence levels
tsSM1_HW_Mensual_forCII <- forecast(tsSM1_HW_MensualII, h=53, level=c(10,25))
## Plot only the forecasted area
plot(tsSM1_HW_Mensual_forCII, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico", main="Predicción para el año 2010", start(2010))

Lavadero
## Seasonal adjusting sub-meter 3 by subtracting the seasonal component & plot
tsSM2_Mensual_AdjustedII <- tsSM2_Mensual ## - Componentes_SM2_Mensual$seasonal
autoplot(tsSM2_Mensual_AdjustedII)

No arece un muelle.

## Volvemos a descomponer la serie
plot(decompose(tsSM2_Mensual_AdjustedII)) 

TENDENCIA CLARAMENTE DECRECIENTE. Tambien observamos esacionalidad

Suavizado exponencial de HoltWinters
## Holt Winters Exponential Smoothing & Plot
tsSM2_HW_MensualII <- HoltWinters(tsSM2_Mensual_AdjustedII, beta=TRUE, gamma=TRUE)
plot(tsSM2_HW_MensualII,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Pronóstico de HoltWinters para el año 2010

Vamos a predecir las próximas 12 observaciones.

## HoltWinters forecast & plot
tsSM2_HW_Mensual_forII <- forecast(tsSM2_HW_MensualII, h=53)
plot(tsSM2_HW_Mensual_forII , ylab= "Vatios-Hora", xlab="Fecha - Lavadero",
     main="Predicción del consumo energético del lavadero \ Año 2010" ) 

## Forecast HoltWinters with diminished confidence levels
tsSM2_HW_Mensual_forCII <- forecast(tsSM2_HW_MensualII, h=53, level=c(10,25))
## Plot only the forecasted area
plot(tsSM2_HW_Mensual_forCII, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico", main="Predicción para el año 2010", start(2010))

AC y termo eléctrico
## Seasonal adjusting sub-meter 3 by subtracting the seasonal component & plot
tsSM3_Mensual_AdjustedII <- tsSM3_Mensual # - Componentes_SM3_Mensual$seasonal
autoplot(tsSM3_Mensual_AdjustedII)

No parece un muelle.

## Volvemos a descomponer la serie
plot(decompose(tsSM3_Mensual_AdjustedII))

Suavizado exponencial de HoltWinters
## Holt Winters Exponential Smoothing & Plot
tsSM3_HW_MensualII <- HoltWinters(tsSM3_Mensual_AdjustedII, beta=TRUE, gamma=TRUE)
plot(tsSM3_HW_MensualII,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Pronóstico de HoltWinters para el año 2010

Vamos a predecir las próximas 12 observaciones.

## HoltWinters forecast & plot
tsSM3_HW_Mensual_forII <- forecast(tsSM3_HW_MensualII, h=53)
plot(tsSM3_HW_Mensual_forII , ylab= "Vatios-Hora", xlab="Fecha - Lavadero",
     main="Predicción del consumo energético del AC y termo eléctrico \ Año 2010" ) 

## Forecast HoltWinters with diminished confidence levels
tsSM3_HW_Mensual_forCII <- forecast(tsSM3_HW_MensualII, h=53, level=c(10,25))
## Plot only the forecasted area
plot(tsSM3_HW_Mensual_forCII, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico", main="Predicción para el año 2010", start(2010))

Comparación de los coeficientes y resultados de cada predicción. Datos con estacionalidad
RMSE_SM1_GranMensualForII<-accuracy(tsSM1_HW_Mensual_forCII   ,DatosRealesSeriesMensuales$Sub_metering_1)[4]
RMSE_SM2_GranMensualForII<-accuracy(tsSM2_HW_Mensual_forCII   ,DatosRealesSeriesMensuales$Sub_metering_2)[4]
RMSE_SM3_GranMensualForII<-accuracy(tsSM3_HW_Mensual_forCII   ,DatosRealesSeriesMensuales$Sub_metering_3)[4]


ResEvolMensual2SinTendenciaConEstacionalidad<-
cbind.data.frame(
      Submedidor=c("Cocina","Lavadero","AC y TermoE"),
      RMSE = c(RMSE_SM1_GranMensualForII,
               RMSE_SM2_GranMensualForII,
               RMSE_SM3_GranMensualForII)
   
    
)

ResEvolMensual2SinTendenciaConEstacionalidad %>% kable(booktabs=TRUE) %>% kable_styling(latex_options = "striped")
Submedidor RMSE
Cocina 16362.16
Lavadero 27223.90
AC y TermoE 165249.73

Comparación de resultados. Predicción antes y despues de eliminar estacionalidad

Comparacion<-cbind.data.frame(
  ResEvolMensual2_SinEstacionalidadConTendencia,
  ResEvolMensual2SinTendenciaConEstacionalidad[,2])

colnames(Comparacion) <- c("Submedidor","RMSE Con Estacionalidad","RMSE Sin estacionalidad")
Comparacion %>% kable(booktabs=TRUE) %>% kable_styling(latex_options = "striped")
Submedidor RMSE Con Estacionalidad RMSE Sin estacionalidad
Cocina 45122.42 16362.16
Lavadero 44294.27 27223.90
AC y TermoE 318158.49 165249.73

Resultados más fiables y con menor error trás haber eliminado la estacionalidad de la serie. También así podemos ver la tendencia.

#p1<-autoplot(tsSM1_Mensual_Adjusted)
#p2<-autoplot(tsSM1_Mensual_AdjustedII)
#grid.arrange(
#             grobs = list(p1,p2),
#           nrow = 2,
#             top="Sin estacionalidad",bottom="Con estacionalidad"
#             )

Evolución diaria del consumo por meses. Mes de Enero

Serie temporal y visualización

En este caso, crearemos una serie para cada mes, ya que la frecuencia de los datos no será la misma, al tener cada mes un número de días diferente.

par(mfrow=c(2,2))
seriesDiariasEnero <- Granularidad_dias %>% filter(year != 2006 & year != 2010 & month == 1)
## Creamos un objeto TS
   tsSM1_DiariaEnero<- ts(seriesDiariasEnero$Sub_metering_1, frequency=31, start=c(2007,1))
   
   tsSM2_DiariaEnero<- ts(seriesDiariasEnero$Sub_metering_2, frequency=31, start=c(2007,1))
   
   tsSM3_DiariaEnero<- ts(seriesDiariasEnero$Sub_metering_3, frequency=31, start=c(2007,1))
   
tsSM_GAP_DiariaEnero<- ts(seriesDiariasEnero$Global_active_power, frequency=31, start=c(2007,1))
## Plot sub-meter 3 with autoplot - add labels, color
autoplot(tsSM1_DiariaEnero, ts.colour = 'red', xlab = "Día", ylab = "Energía (Varios-hora)", main = "Evolución diaria del consumo de energía de la cocina en el mes de Enero \n Años 2007-2009")

autoplot(tsSM2_DiariaEnero, xlab = "Día", ylab = "Energía (Varios-hora)", main = "Evolución diaria del consumo de energía de la lavandería en el mes de Enero \n Años 2007-2009")

autoplot(tsSM3_DiariaEnero, ts.colour = 'blue', xlab = "Día", ylab = "Energía (Varios-hora)", main = "Evolución diaria del consumo de energía del aire acondicionado y termo eléctrico en el mes de Enero \n Años 2007-2009")

autoplot(tsSM_GAP_DiariaEnero, ts.colour = 'green', xlab = "Día", ylab = "Energía (Varios-hora)", main = "Evolución diaria del consumo de energía en el mes de Enero \n Años 2007-2009")

Parece que no existen muchos patrones que nos podrían ayudar a hacer predicciones, no serían muy fiables. (gráfica muelle)

Predicción del año 2010 antes de descomponer la serie

Vamos a plicar regresión lineal a cada serie temporal, y veremos el valor de \(R^2\) y del error cuadrático medio.

SM1. Cocina

Modelo

fitSM1_DiariaEnero <- tslm(tsSM1_DiariaEnero ~ trend + season) 
summary(fitSM1_DiariaEnero)

Call:
tslm(formula = tsSM1_DiariaEnero ~ trend + season)

Residuals:
    Min      1Q  Median      3Q     Max 
-3996.9 -1161.0  -265.1   822.5  6806.5 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)  771.999   1295.981   0.596   0.5536  
trend          9.479      8.968   1.057   0.2947  
season2     -279.146   1787.319  -0.156   0.8764  
season3     1429.042   1787.386   0.800   0.4271  
season4      145.896   1787.499   0.082   0.9352  
season5      501.083   1787.656   0.280   0.7802  
season6      -22.063   1787.859  -0.012   0.9902  
season7      855.125   1788.106   0.478   0.6342  
season8      514.979   1788.399   0.288   0.7744  
season9      399.500   1788.736   0.223   0.8240  
season10     474.687   1789.118   0.265   0.7917  
season11     865.541   1789.545   0.484   0.6304  
season12    1036.062   1790.017   0.579   0.5649  
season13    3653.916   1790.533   2.041   0.0456 *
season14    2432.437   1791.095   1.358   0.1794  
season15     330.291   1791.701   0.184   0.8544  
 [ reached getOption("max.print") -- omitted 16 rows ]
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2189 on 61 degrees of freedom
Multiple R-squared:  0.2234,    Adjusted R-squared:  -0.1713 
F-statistic: 0.566 on 31 and 61 DF,  p-value: 0.957

Resultados:

  • Valor del coeficiente de determinación: 0.2234. Es un valor muy bajo, el modelo únicamente explica el 22% de la variabilidad total.

  • Hay un coeficiente significativamente no nulo.

Predicción IC

Predicción para el mes de Enero de 2010

forecastfitSM1_DiariaEnero <- forecast(fitSM1_DiariaEnero, h=31, level=c(80,90))
## h=12 para predecir el año 2010
## Plot sub-meter 3 forecast, limit y and add labels
plot(forecastfitSM1_DiariaEnero, ylab= "Watt-Hours", xlab="Fecha",
     main="Predicción del consumo energético de Enero del año 2010 \n  a través de un modelo de regresión")

SM2. Lavadero

Modelo

fitSM2_DiariaEnero <- tslm(tsSM2_DiariaEnero ~ trend + season) 
summary(fitSM2_DiariaEnero)

Call:
tslm(formula = tsSM2_DiariaEnero ~ trend + season)

Residuals:
    Min      1Q  Median      3Q     Max 
-4045.7 -1317.7  -282.9  1418.7  4220.7 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)  
(Intercept)  2217.675   1483.490   1.495   0.1401  
trend          -3.261     10.265  -0.318   0.7518  
season2       696.594   2045.917   0.340   0.7347  
season3     -1074.479   2045.995  -0.525   0.6014  
season4      3980.449   2046.123   1.945   0.0563 .
season5     -1766.291   2046.304  -0.863   0.3914  
season6     -1070.697   2046.535  -0.523   0.6027  
season7      1986.897   2046.819   0.971   0.3355  
season8      -674.175   2047.153  -0.329   0.7430  
season9      -885.581   2047.539  -0.433   0.6669  
season10      900.346   2047.977   0.440   0.6618  
season11      927.940   2048.465   0.453   0.6522  
season12     -943.799   2049.005  -0.461   0.6467  
season13     -847.205   2049.597  -0.413   0.6808  
season14     -450.611   2050.239  -0.220   0.8268  
season15      594.649   2050.933   0.290   0.7728  
 [ reached getOption("max.print") -- omitted 16 rows ]
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2506 on 61 degrees of freedom
Multiple R-squared:  0.2859,    Adjusted R-squared:  -0.07693 
F-statistic: 0.788 on 31 and 61 DF,  p-value: 0.7627

Resultados:

  • Valor del coeficiente de determinación: 0.2859. Es un valor bajo, el modelo únicamente explica el 29% de la variabilidad total.

  • Hay tres coeficientes significativamente no nulos.

Predicción IC

Predicción para los próximos 12 meses (año 2010)

forecastfitSM2_DiariaEnero <- forecast(fitSM2_DiariaEnero, h=31, level=c(80,90))
## h=12 para predecir el año 2010
## Plot sub-meter 3 forecast, limit y and add labels
plot(forecastfitSM2_DiariaEnero, ylab= "Watt-Hours", xlab="Fecha",
     main="Predicción del consumo energético del año 2010 \n  a través de un modelo de regresión")

SM3. AC y termo

Modelo

fitSM3_DiariaEnero <- tslm(tsSM3_DiariaEnero~ trend + season) 
summary(fitSM3_DiariaEnero)

Call:
tslm(formula = tsSM3_DiariaEnero ~ trend + season)

Residuals:
    Min      1Q  Median      3Q     Max 
-5223.2 -1643.5  -128.2  1968.5  6090.7 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 7359.3708  1811.1381   4.063  0.00014 ***
trend         -0.8762    12.5325  -0.070  0.94449    
season2     1941.2095  2497.7856   0.777  0.44006    
season3     -337.9143  2497.8799  -0.135  0.89284    
season4     1701.2952  2498.0371   0.681  0.49842    
season5     2241.8380  2498.2572   0.897  0.37305    
season6     1754.3809  2498.5401   0.702  0.48525    
season7     3264.9237  2498.8858   1.307  0.19627    
season8     5867.7999  2499.2943   2.348  0.02215 *  
season9     2021.6760  2499.7656   0.809  0.42180    
season10    3640.5522  2500.2996   1.456  0.15051    
season11    1613.7617  2500.8963   0.645  0.52117    
season12    2887.3045  2501.5556   1.154  0.25292    
season13    4624.5140  2502.2776   1.848  0.06943 .  
season14    3936.0569  2503.0621   1.572  0.12101    
season15    5955.9331  2503.9090   2.379  0.02052 *  
 [ reached getOption("max.print") -- omitted 16 rows ]
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3059 on 61 degrees of freedom
Multiple R-squared:  0.3269,    Adjusted R-squared:  -0.01514 
F-statistic: 0.9557 on 31 and 61 DF,  p-value: 0.5438

Resultados:

  • Valor del coeficiente de determinación:0.3269. Es un valor que es muy bajo, el modelo únicamente explica el 32% de la variabilidad total.

  • Hay tres coeficientes significativamente no nulos.

Predicción IC

Predicción para los próximos 12 meses (año 2010)

forecastfitSM3_DiariaEnero <- forecast(fitSM3_DiariaEnero, h=31, level=c(80,90))
## h=12 para predecir el año 2010
## Plot sub-meter 3 forecast, limit y and add labels
plot(forecastfitSM3_DiariaEnero, ylab= "Watt-Hours", xlab="Fecha",
     main="Predicción del consumo energético del mes de Enero del año 2010 \n  a través de un modelo de regresión")

Se prevee una bajada del consumo

Comparación de los coeficientes y resultados de cada modelo

DatosRealesSeriesEnero<- Granularidad_meses %>% filter(year == 2010 )

RMSE_SM1_GranMensualEnero<-accuracy(forecastfitSM1_DiariaEnero  ,DatosRealesSeriesEnero$Sub_metering_1)[4]
R2_SM1_GranMensualEnero<-0.5796


RMSE_SM2_GranMensualEnero<-accuracy(forecastfitSM2_DiariaEnero   ,DatosRealesSeriesEnero$Sub_metering_2)[4]
R2_SM2_GranMensualEnero<-0.6541



RMSE_SM3_GranMensualEnero<-
  accuracy(
    forecastfitSM3_DiariaEnero ,DatosRealesSeriesEnero$Sub_metering_3)[4]
R2_SM3_GranMensualEnero<-0.7577

ResEvolMensualEnero<-
cbind.data.frame(
      Submedidor=c("Cocina","Lavadero","AC y TermoE"),
      RMSE = c(RMSE_SM1_GranMensualEnero,
               RMSE_SM2_GranMensualEnero,
               RMSE_SM3_GranMensualEnero),
     # MASE = c(MASE_SM1_GranMensual,MASE_SM2_GranMensual,MASE_SM3_GranMensual),
      R2 = c(R2_SM1_GranMensualEnero,
             R2_SM2_GranMensualEnero,
             R2_SM3_GranMensualEnero)
)

ResEvolMensualEnero %>% kable(booktabs=TRUE) %>% kable_styling(latex_options = "striped")
Submedidor RMSE R2
Cocina 42403.33 0.5796
Lavadero 46190.44 0.6541
AC y TermoE 310231.12 0.7577
# MASE: Mean absolute error
# RMSE: root mean scuare error. Diferencia entre los valores predichos y los observados.

El RMSE del modelo con mejor coeficiente de determinación es el más alto. Pero no hay que olvidar que el consumo eléctrico para el submedidor correspondiente al AC y termo eléctrico es el que más energía consume de todos con una diferencia muy grande.

Forecasting descomponiendo la serie (así ok)

Descomponer: quitar tendencia

Descomposición de la serie y visualización

Submetering 1

Se aprecia componente estacional: oscilacione que se repiten cada año. La componente tendencia tambien se aprecia. Tendencia decreciente.

 tsSM1_DiariaEnero %>%
  stl( #t.window=13,
      s.window="periodic", 
      robust=TRUE)  %>% summary()
 Call:
 stl(x = ., s.window = "periodic", robust = TRUE)

 Time.series components:
    seasonal             trend             remainder        
 Min.   :-1324.873   Min.   : 871.8631   Min.   :-6121.092  
 1st Qu.: -761.466   1st Qu.:1309.3572   1st Qu.: -491.953  
 Median : -158.424   Median :1491.6212   Median :  220.447  
 Mean   :    0.000   Mean   :1522.7948   Mean   :  550.818  
 3rd Qu.:  337.611   3rd Qu.:1775.0069   3rd Qu.:  732.488  
 Max.   : 5397.352   Max.   :2077.1908   Max.   : 9934.040  
 IQR:
     STL.seasonal STL.trend STL.remainder data  
     1099.1        465.6    1224.4        1692.0
   %  65.0         27.5      72.4         100.0 

 Weights:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.0000  0.8124  0.9452  0.7973  0.9880  1.0000 

 Other components: List of 5
 $ win  : Named num [1:3] 931 47 31
 $ deg  : Named int [1:3] 0 1 1
 $ jump : Named num [1:3] 94 5 4
 $ inner: int 1
 $ outer: int 15

Remainder no es un muelle. Observamos tendencia creciente y fuerte componente estacional

Componentes_SM1_DiariaEneroSTL<-
  tsSM1_DiariaEnero %>%stl( #t.window=13,
  s.window="periodic",  robust=TRUE)
autoplot(Componentes_SM1_DiariaEneroSTL, ts.colour = 'blue')

plot.ts(tsSM1_Mensual, col = 'gray')
lines(Componentes_SM1_DiariaEneroSTL$time.series[,2], 
      col = "red", lwd = 1, lty = 2)

Submetering 2

 tsSM2_DiariaEnero %>%
  stl( #t.window=13,
      s.window="periodic", 
      robust=TRUE) %>% summary()
 Call:
 stl(x = ., s.window = "periodic", robust = TRUE)

 Time.series components:
    seasonal             trend            remainder        
 Min.   :-1788.921   Min.   :1491.391   Min.   :-5599.893  
 1st Qu.:-1294.941   1st Qu.:1617.887   1st Qu.: -685.829  
 Median : -782.743   Median :1866.286   Median : -122.576  
 Mean   :    0.000   Mean   :1939.768   Mean   :  373.264  
 3rd Qu.:  718.901   3rd Qu.:2339.704   3rd Qu.: 1313.902  
 Max.   : 4203.925   Max.   :2437.168   Max.   : 6567.217  
 IQR:
     STL.seasonal STL.trend STL.remainder data  
     2013.8        721.8    1999.7        3494.0
   %  57.6         20.7      57.2         100.0 

 Weights:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.0000  0.8067  0.9452  0.7789  0.9886  1.0000 

 Other components: List of 5
 $ win  : Named num [1:3] 931 47 31
 $ deg  : Named int [1:3] 0 1 1
 $ jump : Named num [1:3] 94 5 4
 $ inner: int 1
 $ outer: int 15
Componentes_SM2_DiariaEneroSTL<- tsSM2_DiariaEnero %>%stl( #t.window=13,
  s.window="periodic",  robust=TRUE)
autoplot(Componentes_SM2_DiariaEneroSTL, ts.colour = 'blue')

plot.ts(tsSM2_DiariaEnero, col = 'gray')
lines(Componentes_SM2_DiariaEneroSTL$time.series[,2], 
      col = "red", lwd = 1, lty = 2)

Submetering 3

## Descomponer la serie
Componentes_SM3_DiariaEnero <- decompose(tsSM3_DiariaEnero)
## Plot
plot(Componentes_SM3_DiariaEnero)

## Check summary statistics for decomposed sub-meter 2
summary(Componentes_SM3_DiariaEnero)
         Length Class  Mode     
x        93     ts     numeric  
seasonal 93     ts     numeric  
trend    93     ts     numeric  
random   93     ts     numeric  
figure   31     -none- numeric  
type      1     -none- character

Se aprecia componente estacional: oscilacione que se repiten cada año. La componente tendencia tambien se aprecia. Tendencia creciente

tsSM3_DiariaEnero %>%
  stl( #t.window=13,
      s.window="periodic", 
      robust=TRUE)  %>% summary()
 Call:
 stl(x = ., s.window = "periodic", robust = TRUE)

 Time.series components:
    seasonal             trend             remainder        
 Min.   :-5596.817   Min.   : 9448.476   Min.   :-8001.466  
 1st Qu.:-1429.752   1st Qu.:10055.012   1st Qu.:-1298.672  
 Median :  135.662   Median :10467.966   Median : -172.626  
 Mean   :    0.000   Mean   :10408.787   Mean   :   31.439  
 3rd Qu.: 1423.036   3rd Qu.:10742.199   3rd Qu.: 1502.564  
 Max.   : 3704.229   Max.   :11008.213   Max.   : 8713.785  
 IQR:
     STL.seasonal STL.trend STL.remainder data  
     2852.8        687.2    2801.2        3550.0
   %  80.4         19.4      78.9         100.0 

 Weights:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.0000  0.8394  0.9452  0.8553  0.9931  1.0000 

 Other components: List of 5
 $ win  : Named num [1:3] 931 47 31
 $ deg  : Named int [1:3] 0 1 1
 $ jump : Named num [1:3] 94 5 4
 $ inner: int 1
 $ outer: int 15
Componentes_SM3_DiariaEneroSTL<- tsSM3_DiariaEnero %>%stl( #t.window=13,
  s.window="periodic",  robust=TRUE)
autoplot(Componentes_SM3_DiariaEneroSTL, ts.colour = 'blue')

plot.ts(tsSM3_DiariaEnero, col = 'gray')
lines(Componentes_SM3_DiariaEnero$time.series[,2], 
      col = "red", lwd = 1, lty = 2)

Predicción con Holt-Winters (suavizado exponencial) para el año 2010

Sin estacionalidad

Sería interesante verlo también con estacionalidad, para ver si se mantiene o no la tendencia.

Tendencia de predicción (crecimiento o decrecimiento del consumo) y predicción teniendo en cuenta las componentes

Evolución mensual del consumo (sin estacionalidad). Años 2007-2009

Cocina
tsSM1_DiariaEnero_detrended <- tsSM1_DiariaEnero - Componentes_SM1_DiariaEneroSTL$time.series[,2]
plot.ts(tsSM1_DiariaEnero_detrended, 
        ylab = "Energía (Vatios-hora)", 
        main = 'Seasonal variability of energy comsumption',
        cex.main = 0.85)

par(mfrow = c(1,2))
  plot(Componentes_SM1_DiariaEneroSTL$time.series[,3], col = "blue", main = 'Remainder', ylab = "")
qqnorm(Componentes_SM1_DiariaEneroSTL$time.series[,3])
qqline(Componentes_SM1_DiariaEneroSTL$time.series[,3])

shapiro.test(Componentes_SM1_DiariaEneroSTL$time.series[,3])

    Shapiro-Wilk normality test

data:  Componentes_SM1_DiariaEneroSTL$time.series[, 3]
W = 0.78722, p-value = 2.838e-10

Los resíduos no son solo ruido. Aún hay estacionalidad

Volvemos a descomponer la serie desestacionalizada una vez.

## Volvemos a descomponer la serie
plot(stl(tsSM1_DiariaEnero_detrended ,s.window = "periodic",robust = TRUE))

La componente estacional se mueve en un ranfo de \(10^{-11}\), podemos asumir que se ha eliminado esta componente.

Suavizado exponencial de HoltWinters
## Holt Winters Exponential Smoothing & Plot
tsSM1_HW_DiariaEnero<- HoltWinters(tsSM1_DiariaEnero_detrended ,
                              #  tsSM1_Mensual_Adjusted, 
  beta=TRUE, gamma=TRUE)
# optimiza
plot(tsSM1_HW_DiariaEnero,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Pronóstico de HoltWinters para el año 2010

Vamos a predecir las próximas 12 observaciones.

## HoltWinters forecast & plot
tsSM1_HW_DiariaEnero_for <- forecast(tsSM1_HW_DiariaEnero, h=31)
plot(tsSM1_HW_DiariaEnero_for , ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico",
     main="Predicción del consumo energético de la cocina \ Año 2010" ) 

## Forecast HoltWinters with diminished confidence levels
tsSM1_HW_DiariaEnero_forC <- forecast(tsSM1_HW_DiariaEnero, h=31, level=c(10,25))
## Plot only the forecasted area
plot(tsSM1_HW_DiariaEnero_forC, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico", main="Predicción para el año 2010", start(2010))

Lavadero
tsSM2_DiariaEnero_detrended <- tsSM2_DiariaEnero - Componentes_SM2_DiariaEneroSTL$time.series[,2]
plot.ts(tsSM2_DiariaEnero_detrended, 
        ylab = "Energía (Vatios-hora)", 
        main = 'Seasonal variability of energy comsumption',
        cex.main = 0.85)

par(mfrow = c(1,2))
  plot(Componentes_SM2_DiariaEneroSTL$time.series[,3], col = "blue", main = 'Remainder', ylab = "")
qqnorm(Componentes_SM2_DiariaEneroSTL$time.series[,3])
qqline(Componentes_SM2_DiariaEneroSTL$time.series[,3])

shapiro.test(Componentes_SM2_DiariaEneroSTL$time.series[,3])

    Shapiro-Wilk normality test

data:  Componentes_SM2_DiariaEneroSTL$time.series[, 3]
W = 0.88372, p-value = 5.75e-07

Los resíduos no siguen una distribución normal. La descomposición no es buena.

## Volvemos a descomponer la serie
plot(stl(tsSM2_DiariaEnero_detrended ,s.window = "periodic",robust = TRUE))

La componente estacional se mueve en un ranfo de \(10^{-11}\), podemos asumir que se ha eliminado esta componente.

Suavizado exponencial de HoltWinters
## Holt Winters Exponential Smoothing & Plot
tsSM2_HW_DiariaEnero <- HoltWinters(tsSM2_DiariaEnero_detrended ,
                                # tsSM2_Mensual_Adjusted,
                                beta=TRUE, gamma=TRUE)
plot(tsSM2_HW_DiariaEnero,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Pronóstico de HoltWinters para el año 2010

Vamos a predecir las próximas 12 observaciones.

## HoltWinters forecast & plot
tsSM2_HW_DiariaEnero_for <- forecast(tsSM2_HW_DiariaEnero, h=31)
plot(tsSM2_HW_DiariaEnero_for , ylab= "Vatios-Hora", xlab="Fecha - Lavadero",
     main="Predicción del consumo energético del lavadero \ Año 2010" ) 

## Forecast HoltWinters with diminished confidence levels
tsSM2_HW_DiariaEnero_forC <- forecast(tsSM2_HW_DiariaEnero, h=31, level=c(10,25))
## Plot only the forecasted area
plot(tsSM2_HW_DiariaEnero_forC, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico", main="Predicción para el año 2010", start(2010))

AC y termo eléctrico
tsSM3_DiariaEnero_detrended <- tsSM3_DiariaEnero - Componentes_SM3_DiariaEneroSTL$time.series[,2]
plot.ts(tsSM3_DiariaEnero_detrended, 
        ylab = "Energía (Vatios-hora)", 
        main = 'Seasonal variability of energy comsumption',
        cex.main = 0.85)

par(mfrow = c(1,2))
  plot(Componentes_SM3_DiariaEneroSTL$time.series[,3], col = "blue", main = 'Remainder', ylab = "")
qqnorm(Componentes_SM3_DiariaEneroSTL$time.series[,3])
qqline(Componentes_SM3_DiariaEneroSTL$time.series[,3])

shapiro.test(Componentes_SM3_DiariaEneroSTL$time.series[,3])

    Shapiro-Wilk normality test

data:  Componentes_SM3_DiariaEneroSTL$time.series[, 3]
W = 0.95542, p-value = 0.003018

Los resíduos no siguen una distribución normal. La descomposición no es buena.

## Volvemos a descomponer la serie
plot(stl(tsSM3_DiariaEnero_detrended ,s.window = "periodic",robust = TRUE))

La componente estacional se mueve en un ranfo de \(10^{-11}\), podemos asumir que se ha eliminado esta componente.

Suavizado exponencial de HoltWinters
## Holt Winters Exponential Smoothing & Plot
tsSM3_HW_DiariaEnero <- HoltWinters(tsSM3_DiariaEnero_detrended ,
                              # tsSM3_Mensual_Adjusted,
                                beta=TRUE, gamma=TRUE)
plot(tsSM3_HW_DiariaEnero,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Pronóstico de HoltWinters para el año 2010

Vamos a predecir las próximas 12 observaciones.

## HoltWinters forecast & plot
tsSM3_HW_DiariaEnero_for <- forecast(tsSM3_HW_DiariaEnero, h=31)
plot(tsSM3_HW_DiariaEnero_for , ylab= "Vatios-Hora", xlab="Fecha - Lavadero",
     main="Predicción del consumo energético del AC y termo eléctrico \ Año 2010" ) 

## Forecast HoltWinters with diminished confidence levels
tsSM3_HW_DiariaEnero_forC <- forecast( object = tsSM3_HW_DiariaEnero_for, h=31)
## Plot only the forecasted area
plot(tsSM3_HW_DiariaEnero_forC, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico", main="Predicción para el año 2010", start(2010))

Comparación de los coeficientes y resultados de cada predicción. Datos sin estacionalidad

En las gráficas de predicción, hemos observado que el ajuste no mantiene la tendencia de consumo energético de la serie.

# RMSE_SM1_GranDiariaEneroFor<-accuracy(tsSM1_HW_DiariaEnero_forC   # ,seriesDiariasEnero$Sub_metering_1)[4]
# RMSE_SM2_GranDiariaEneroFor<-accuracy(tsSM2_HW_DiariaEnero_forC   # ,seriesDiariasEnero$Sub_metering_2)[4]
# RMSE_SM3_GranDiariaEneroFor<-accuracy(tsSM3_HW_DiariaEnero_forC   # ,seriesDiariasEnero$Sub_metering_3)[4]
# 
# 
# ResEvolDiariaEnero2_SinEstacionalidadConTendencia <-
# cbind.data.frame(
#       Submedidor=c("Cocina","Lavadero","AC y TermoE"),
#       RMSE = c(RMSE_SM1_GranDiariaEneroFor,
#                RMSE_SM2_GranDiariaEneroFor,
#                RMSE_SM3_GranDiariaEneroFor)
# )
# ResEvolDiariaEnero2_SinEstacionalidadConTendencia %>% kable(booktabs=TRUE) %>% # kable_styling(latex_options = "striped")

Evolución mensual del consumo (con estacionalidad). Años 2007-2009

Vamos a ver una predicción del consumo enegético de cada submedidor con el método de HW. En este caso, no eliminaremos la estacionalidad, para ver así si se mantiene la tendencia de consumo.

Cocina
tsSM1_DiariaEnero_AdjustedII <- tsSM1_DiariaEnero # - Componentes_SM1_Mensual$seasonal
autoplot(tsSM1_DiariaEnero_AdjustedII)

Parece un muelle.

## Volvemos a descomponer la serie
plot(decompose(tsSM1_DiariaEnero_AdjustedII))

La componente estacional se mueve en un rango de \(10^{-11}\), podemos asumir que se ha eliminado esta componente.

Suavizado exponencial de HoltWinters
## Holt Winters Exponential Smoothing & Plot
tsSM1_HW_DiariaEneroII <- HoltWinters(tsSM1_DiariaEnero_AdjustedII, beta=TRUE, gamma=TRUE)
# optimiza
plot(tsSM1_HW_DiariaEneroII,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Pronóstico de HoltWinters para el año 2010

Vamos a predecir las próximas 12 observaciones.

## HoltWinters forecast & plot
tsSM1_HW_DiariaEnero_forII <- forecast(tsSM1_HW_DiariaEneroII, h=31)
plot(tsSM1_HW_DiariaEnero_forII , ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico",
     main="Predicción del consumo energético de la cocina \ Año 2010" ) 

## Forecast HoltWinters with diminished confidence levels
tsSM1_HW_DiariaEnero_forCII <- forecast(tsSM1_HW_DiariaEneroII, h=31, level=c(10,25))
## Plot only the forecasted area
plot(tsSM1_HW_DiariaEnero_forCII, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico", main="Predicción para el año 2010", start(2010))

Lavadero
## Seasonal adjusting sub-meter 3 by subtracting the seasonal component & plot
tsSM2_DiariaEnero_AdjustedII <- tsSM2_DiariaEnero ## - Componentes_SM2_Mensual$seasonal
autoplot(tsSM2_DiariaEnero_AdjustedII)

No arece un muelle.

## Volvemos a descomponer la serie
plot(decompose(tsSM2_DiariaEnero_AdjustedII)) 

TENDENCIA CLARAMENTE DECRECIENTE. Tambien observamos esacionalidad

Suavizado exponencial de HoltWinters
## Holt Winters Exponential Smoothing & Plot
tsSM2_HW_DiariaEneroII <- HoltWinters(tsSM2_DiariaEnero_AdjustedII, beta=TRUE, gamma=TRUE)
plot(tsSM2_HW_DiariaEneroII,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Pronóstico de HoltWinters para el año 2010

Vamos a predecir las próximas 12 observaciones.

## HoltWinters forecast & plot
tsSM2_HW_DiariaEnero_forII <- forecast(tsSM2_HW_DiariaEneroII, h=31)
plot(tsSM2_HW_DiariaEnero_forII , ylab= "Vatios-Hora", xlab="Fecha - Lavadero",
     main="Predicción del consumo energético del lavadero \ Año 2010" ) 

## Forecast HoltWinters with diminished confidence levels
tsSM2_HW_DiariaEnero_forCII <- forecast(tsSM2_HW_DiariaEneroII, h=31, level=c(10,25))
## Plot only the forecasted area
plot(tsSM2_HW_DiariaEnero_forCII, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico", main="Predicción para el año 2010", start(2010))

AC y termo eléctrico
## Seasonal adjusting sub-meter 3 by subtracting the seasonal component & plot
tsSM3_DiariaEnero_AdjustedII <- tsSM3_DiariaEnero # - Componentes_SM3_Mensual$seasonal
autoplot(tsSM3_DiariaEnero_AdjustedII)

No parece un muelle.

## Volvemos a descomponer la serie
plot(decompose(tsSM3_DiariaEnero_AdjustedII))

Suavizado exponencial de HoltWinters
## Holt Winters Exponential Smoothing & Plot
tsSM3_HW_DiariaEneroII <- HoltWinters(tsSM3_DiariaEnero_AdjustedII, beta=TRUE, gamma=TRUE)
plot(tsSM3_HW_DiariaEneroII,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Pronóstico de HoltWinters para el año 2010

Vamos a predecir las próximas 12 observaciones.

## HoltWinters forecast & plot
     tsSM3_HW_DiariaEnero_forII <- forecast(tsSM3_HW_DiariaEneroII, h=31)
plot(tsSM3_HW_DiariaEnero_forII , ylab= "Vatios-Hora", xlab="Fecha - Lavadero",
     main="Predicción del consumo energético del AC y termo eléctrico \ Año 2010" ) 

## Forecast HoltWinters with diminished confidence levels
     tsSM3_HW_DiariaEnero_forCII <- forecast(tsSM3_HW_DiariaEneroII, h=31, level=c(10,25))
plot(tsSM3_HW_DiariaEnero_forCII, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico", main="Predicción para el año 2010", start(2010))

Comparación de los coeficientes y resultados de cada predicción. Datos con estacionalidad
# RMSE_SM1_GranDiariaEneroForII<-accuracy(tsSM1_HW_DiariaEnero_forCII   # ,seriesDiariasEnero$Sub_metering_1)[4]
# RMSE_SM2_GranDiariaEneroForII<-accuracy(tsSM2_HW_DiariaEnero_forCII   # ,seriesDiariasEnero$Sub_metering_2)[4]
# RMSE_SM3_GranDiariaEneroForII<-accuracy(tsSM3_HW_DiariaEnero_forCII   # ,seriesDiariasEnero$Sub_metering_3)[4]
# 
# 
# ResEvolDiariaEnero2SinTendenciaConEstacionalidad<-
# cbind.data.frame(
#       Submedidor=c("Cocina","Lavadero","AC y TermoE"),
#       RMSE = c(RMSE_SM1_GranDiariaEneroForII,
#                RMSE_SM2_GranDiariaEneroForII,
#                RMSE_SM3_GranDiariaEneroForII)
#    
#     
# )
# 
# ResEvolDiariaEnero2SinTendenciaConEstacionalidad %>% kable(booktabs=TRUE) %>% # kable_styling(latex_options = "striped")

Comparación de resultados. Predicción antes y despues de eliminar estacionalidad

# Comparacion<-cbind.data.frame(
#   ResEvolDiariaEnero2_SinEstacionalidadConTendencia,
#   ResEvolDiariaEnero2SinTendenciaConEstacionalidad[,2])
# 
# colnames(Comparacion) <- c("Submedidor","RMSE Con Estacionalidad","RMSE Sin estacionalidad")
# Comparacion %>% kable(booktabs=TRUE) %>% kable_styling(latex_options = "striped")

Resultados más fiables y con menor error trás haber eliminado la estacionalidad de la serie. También así podemos ver la tendencia.

#p1<-autoplot(tsSM1_Mensual_Adjusted)
#p2<-autoplot(tsSM1_Mensual_AdjustedII)
#grid.arrange(
#             grobs = list(p1,p2),
#           nrow = 2,
#             top="Sin estacionalidad",bottom="Con estacionalidad"
#             )
library(dplyr)
a<-as.data.frame(Granularidad_dias %>% select("energia2","Date") %>% filter(Date > "2008-12-16") )

library(plotly)

 plot_ly(x=~a$Date,
            y = ~a[,3],
            type = 'scatter',mode = 'lines',name="Plot",
            line=list(color='rgb(199, 0, 57 )'),marker = list(color =   'rgb(199, 0, 57  )')
    )
library(tidyr)
DatosRango<-Granularidad_dias%>%
  pivot_longer(cols = c(7,8,9,10,11),
               names_to = "Energia",
values_to = "valorEnergetico") %>% 
  select(
  Date , Energia , valorEnergetico
)

DatosRango %>% filter(Energia == "Sub_metering_1", Date > "2008-12-12")
# A tibble: 715 × 5
# Groups:   year, month [24]
    year month Date                Energia        valorEnergetico
   <dbl> <dbl> <dttm>              <chr>                    <dbl>
 1  2008    12 2008-12-12 01:00:00 Sub_metering_1            1157
 2  2008    12 2008-12-13 01:00:00 Sub_metering_1            2186
 3  2008    12 2008-12-14 01:00:00 Sub_metering_1            3742
 4  2008    12 2008-12-15 01:00:00 Sub_metering_1             592
 5  2008    12 2008-12-16 01:00:00 Sub_metering_1             586
 6  2008    12 2008-12-17 01:00:00 Sub_metering_1            3680
 7  2008    12 2008-12-18 01:00:00 Sub_metering_1            1125
 8  2008    12 2008-12-19 01:00:00 Sub_metering_1            2109
 9  2008    12 2008-12-20 01:00:00 Sub_metering_1            1134
10  2008    12 2008-12-21 01:00:00 Sub_metering_1             922
# … with 705 more rows